Genotyping and inflated type I error rate in genome-wide association case/control studies

Calling Genotypes

where ψ(·) is the multivariate normal density. The three mean vectors, $\mu \equiv \{{\overrightarrow{\mu }}_{AA},{\overrightarrow{\mu }}_{AB},{\overrightarrow{\mu }}_{BB}\}$, variance matrices, Σ ≡ {Σ _AA , Σ _AB , Σ _BB }, and probabilities, p ≡ {p _AA , p _AB , p _BB }, correspond to the three possible genotypes, AA, AB, and BB. Define $\theta \equiv \{\overrightarrow{\mu },\Sigma ,\overrightarrow{p}\}$ and later, we let Φ(·) and Φ ^M (·) be the cumulative distribution for a normal variable and a mixture of normals. The third step is to estimate θ. Some algorithms, such as RLMM, use a training data set, where the genotypes are known. Other algorithms, such as SNiPer-HD, use no training data, and find the best estimates that describe their experimental values. The fourth step is to assign a genotype, ${\hat{G}}_{ij}$, to SNP j in subject i. Often, for a given value of S _i , the assigned genotype maximizes a similarity function: ${\hat{G}}_i$ = argmax _g D(g, S _i |θ). The similarity function is usually a modified version of one of the following three quantities: Mahalanobis distance: $-(S_i-{\overrightarrow{\mu }}_g){\Sigma }_g^{-1}{(S_i-{\overrightarrow{\mu }}_g)}^t$, Unweighted Likelihood: $\varphi (S_i|{\overrightarrow{\mu }}_g,{\Sigma }_g)$, or the Weighted Likelihood: p _g $\varphi (S_i|{\overrightarrow{\mu }}_g,{\Sigma }_g)$. The similarity function is also modified to ensure monotonicity of assignment. When we let S _i = {M _i , A _i }, we force D(AB, S _i |θ) = D(BB, S _i |θ) = 0 if M _i ≤ μ _AA , D(BB, S _i |θ) = 0 if μ _AA ≤ M _i ≤ μ _AB , D(AA, S _i |θ) = 0 if μ _AB ≤ M _i ≤ · μ _BB , and D(AB, S _i |θ) = D(AA, S _i |θ) = 0 if M _i ≥ μ _BB , where μ _g , in this case, is the mean of M _i when G _i = g. This modification is standard in calling algorithms. As we do not know the true value of the parameters in experiments, we replace D(g, S _i |θ) by D(g, S _i |$\hat{\theta }$). A subject's SNP may not be called, or assigned a missing value, if the difference or ratio between D(${\hat{G}}_i$, S _i |$\hat{\theta }$) and D(g_2i, S _i |$\hat{\theta }$) is not large enough, where g_2iis the genotype with the second largest value of D(g, S _i |$\hat{\theta }$). A SNP may be omitted from further study if too many values were set to missing. Table 1 describes the details of the four steps for popular methods. For many purposes or to understand the details of the method, especially in handling rare alleles, this table will seem an oversimplification. For our purposes here, it highlights the features of interest.

Tests of Association

The current tests of association start by calling genotypes for a given SNP j in a group of subjects with the disease and in a group of controls. They then compare the resulting proportions, ${\hat{P}}^A\equiv \{{\hat{P}}_{AA}^A,{\hat{P}}_{AB}^A,{\hat{P}}_{BB}^A\}$ and ${\hat{P}}^U\equiv \{{\hat{P}}_{AA}^U,{\hat{P}}_{AB}^U,{\hat{P}}_{BB}^U\}$, from these Affected and Unaffected groups using either a Cochran-Armitage test or logistic regression. Here ${\hat{p}}_g^q\equiv \frac{1}{nq}{\displaystyle {\sum }_{i:Q_{i=q}}1({\hat{G}}_i=g)}$, where the indicator function is defined by 1(x) = 1 if x is true, 0 otherwise, Q _i ∈ {A, U} is the disease status for individual i, and n _q is the number of subjects with disease status q. In this manuscript, any 'p-value' from a genotype-based association test will be calculated using ANOVA on the logistic regression model with Q _i and genotype (unordered grouping) as the dependent and independent variables.

Standard tests tend to err anti-conservatively as we will discuss below. We will propose four alterations that can reduce type I error rate, with only a minimal decrease in power. These are the four differences that separate logiCALL from standard methods. The first is based on the observation, which is discussed later, that the likelihood profile of φ ^M (S _i |θ) will have multiple local maxima near the overall maximum. When estimating θ, the EM algorithm converges to multiple solutions. For many of those solutions, the resulting parameter set, ${\hat{\theta }}_{lm}$, satisfies ${\displaystyle {\prod }_i{\varphi }^M(S_i|\hat{\theta })\le }{\displaystyle {\prod }_i{\varphi }^M(S_i|{\hat{\theta }}_{lm})}+\tau n$. For each parameter set satisfying this inequality, we will make a new group of genotype assignments, $\left\{{\hat{G}}_i^{lm}\right\}$. If more than 10% of such assignments disagree with ${\hat{G}}_i$ (τ = 0.06), we label that subject's call as questionable. We also continue the practice of marking calls with small values of D(${\hat{G}}_i$, S _i |$\hat{\theta }$)/D(g_2i, S _i |$\hat{\theta }$) as questionable. The second alteration is that we do not discard questionable calls, an act which can create false positives. Instead, we assign questionable S _i so the proportions of genotypes in the cases and controls are as similar as possible, which is defined as minimizing ${\displaystyle {\sum }_g|{\hat{p}}_g^U-{\hat{p}}_g^A|}$, with the restriction that the final call for subject i must be either ${\hat{G}}_i$ or g_2i. Here, we let g_2ibe the genotype which is either the runner-up in terms of distance or the most common genotype among the dissenting calls, depending on why the genotype was labeled as questionable. The third alteration, which is already incorporated into other programs is to perform the EM algorithm under the null hypothesis that the genotype proportions in the two populations are identical [3]. The fourth is the use of a weighted Mahalonobis distance, which is defined later. Given these changes, logiCALL then compares the estimated genotype proportions in cases and controls using logistic regression. Note that none of the changes affect calls for the vast majority of SNPs.

Clearly, the restricted parameter space, $\{{\Omega }_R:p_{AA}^A,p_{AA}^U,p_{AB}^A=p_{AB}^U,p_{BB}^A=p_{BB}^U\}$ is a subset of the unrestricted parameter space, Ω. In an ideal scenario, the distribution of 2log(LR) would converge to a chi-squared distribution, $F_{{\chi }_2^2}$, with 2 degrees of freedom. Therefore, the 'p-value' from a likelihood ratio-based test will be calculated as 1 - $F_{{\chi }_2^2}$(2log(LR)).

Data Source

To demonstrate the problems of genotyping and compare the genotype- and likelihood ratio-based tests of association, we use three types of data. First, for discussion, we may assume a hypothetical study measuring a one dimensional summary statistic, S _i , for a SNP j with only two possible genotypes, G _i ∈ {0, 1}. Furthermore, to show problems can exist even under the best conditions, where model and truth coincide, we assume that S _i follows a normal distribution given Q _i and G _ij , and that the full distribution can be described by ${\varphi }^M(S_i|\theta )=p_0\varphi (S_i|{\mu }_0,{\sigma }_0^2)+p_1\varphi (S_i|{\mu }_1,{\sigma }_1^2)$.

We compare our two new tests of association to a standard method using simulated data. The standard method mimics the general Bead-Studio approach by a) fitting parameters with the EM algorithm; b) calling genotypes based on the Mahalanobis Distance; c) removing all calls where D(${\hat{G}}_i$, S _i |$\hat{\theta }$)/D(g₂, S _i |$\hat{\theta }$) > 0.5; and d) comparing the two sets of resulting estimates, $\{{\hat{P}}_{AA}^A,{\hat{P}}_{AB}^A,{\hat{P}}_{BB}^A\}$ and $\{{\hat{P}}_{AA}^U,{\hat{P}}_{AB}^U,{\hat{P}}_{BB}^U\}$. We generated 10 simulated datasets, containing 1000 subjects (500 cases, 500 controls) and 303,100 SNPs for each of 18 scenarios. For each gene j and each subject i, we generated a 2D summary statistic (M _ji , A _ji ). The distribution of M _ji depended on genotype. If G _ji = AA, then M _ji ~ 2X - 1, and if G _ji = BB, then M _ji ~1 - 2X, where X ~ beta(α = 3, β = 30). If G _ji = AB, then M _ji also followed a beta distribution, but the parameters varied by SNP, disease status, and scenario. For all SNPs and all subjects, A _ji ~ N (10, 1.5). We generated three types of SNPs, background, shifted, and influential. First, 300,000 background SNPs were generated and included in all 10 × 18 = 180 data sets. For each SNP, a single minor allele frequency was generated from a uniform(0.2, 0.4) distribution and genotype probabilities $({\overrightarrow{p}}^U={\overrightarrow{p}}^A)$ were generated assuming Hardy-Weinberg Equilibrium. Here, E [M _ji |G _ji = AB] = 0 and was independent of disease status. These SNPs, which formed three distinct clusters, can be easily identified and represent a well-behaved group. For 3,000 shifted SNPs, MAF ~uniform(0.2, 0.4) and genotype probabilities $({\overrightarrow{p}}^U={\overrightarrow{p}}^A)$ were generated assuming Hardy-Weinberg Equilibrium. Here, E [M _ji |G _ji = AB] ∈ {-0.739 + 0.2, - 0.739 + 0.3 - 0.739 + 0.5}, where we note E [M _ji |G _ji = AA] = -0.739 and E [M _ji |G _ji = AB, Q _i = A] - E [M _ji |G _ji = AB, Q _i = U] ∈ {0, 0.2}. This group represents difficult to call SNPs. For 100 influential SNPs, MAF ~uniform(0.2, 0.4) and genotype probabilities $({\overrightarrow{p}}^U\ne {\overrightarrow{p}}^A)$ were chosen so that, under a disease prevalence of 0.01 and a model of additive effects, the genotype relative risk for subjects homogeneous for the minor allele, P(Q _i = A|BB)/P (Q _i = A|AA) ∈ {1.5, 2.0, 2.5}. Combing the degree of shift for poor quality SNPs and the effect size of truly associated SNPs, we have a total of 18 scenarios used in our simulation.

The next set of data is from a recent GWAS of Inflammatory Bowel Disease (IBD) that compared 983 subjects with IBD to 1004 subjects without the disease. Using Illumina microarrays, 308,330 SNPs on the autosomal chromosomes were tested. Jewish and non-Jewish cohorts, approximately equal in size, were analyzed separately, a practice continued here. Details have been previously published [17, 18]. Because the overwhelming majority of the SNPs are easy to genotype, as any of the summary statistics neatly divide into three clusters, we chose a 3137 difficult SNP subset where at least two clusters overlap (definition below). Because association was tested separately in Jewish and non-Jewish cohorts, a total of 3137 × 2 = 6274 tests were possible. To demonstrate called genotype dependency, we bootstrapped 40 samples, ignoring case/control status, of 500 subjects for each SNP. A sample would be discarded, and replaced by another random selection, if one or both of the groups were lacking an AA (S _i <-0.7) or a BB (S _i > 0.7) genotype.

We defined difficult SNPs as follows. For SNP j, we first estimated the density, $\hat{f}$(M_j.) nonparametrically using the R function 'density(adjust = 0.3)'. In theory, $\hat{f}$(M_j.) is a mixture of three normal distributions, corresponding to the three genotypes. If the SNP is well-behaved, then the three underlying densities will not overlap, and the empirical density $\hat{f}$(M_j.) will attain minima near 0 in the valleys between μ _jAA and μ _jAB and between μ _jAB and μ _jBB . If either of these minima exceeded 0.2, then at least two of underlying densities overlapped, and that SNP was defined as difficult. To speed the process, we found that approximating the center of the peaks (i.e. μ _jAA , μ _jAB , and μ _jBB ) by the median values of M_j.in the three windows, {M_j.≤ -0.6; -0.3 ≤ M_j.≤ 0.3, M_j.≥ 0.5}, worked well.

Two Genotype Example: Parameters

We choose to use the hypothetical, two-genotype, study, to highlight that the estimated parameters can be inconsistent even in the simplest scenario, where the summary statistic is distributed normally and our fitted model is correct. When dealing with only a single population, we define the parameter p₀ (p₁) to be the probability that a subject's true genotype is 0 (1).

p₀ ≡ P (G _i = 0) and p₁ ≡ P(G _i = 1) (4)

Probabilities of called genotypes implicitly depend on the number of subjects in the sample. We then define the parameter c* to be the point which is equidistant to the two genotype groups when θ is known.

Therefore, μ₀ ≤ c* ≤ μ₁ is defined to be a solution to equation 6

D(G _i = 0, S _i = c*|θ) = D(G _i = 1, S _i = c*|θ) (6)

Therefore, by convergence of the MLE, we know that

lim_{n → ∞}C → _p c* (12)

which will useful for the next section.

Two Genotype Example: $p_0^{*}\ne p_0$

Clearly, when p₀ = p₁ = 0.5, $p_0^{*}$ for all values of $P_M^{*}$. However, when p₀ ≠ p₁, the bias is a non-zero function $P_M^{*}$, and therefore depends on the parameter group, $\{{\mu }_1-{\mu }_0,{\sigma }_1^2/{\sigma }_0^2\}$ (Figure 1). As shown by Figure 1, the bias can be quite large when either $P_M^{*}$ and/or |p₀ - 0.5 | is large.

Next, assume that genotype assignments are based on a likelihood, or weighted likelihood measure. For a value ω ∈ c(0, 1), the probability that φ _g (S _i ) exceeds ω, or $2{\Phi }_g({\varphi }_g^{-1}(\omega ))-1$, changes with ${\sigma }_g^2$, where ${\varphi }_g^{-1}(\omega )$ returns a value greater than μ _g . Therefore, φ₀(s) = φ₁(s) does not imply anything about the relationship between Φ₀(s) = Φ₁(s). We illustrate the potential for bias by a simple example where μ₀ = 0, ${\sigma }_0^2$ = 1, and μ₁ = 1. Let p₀ = p₁ = 0.5. Then, for any value of ${\sigma }_1^2\ne 1,1-\Phi ({\varphi }_0^{-1}(\omega ))\ne \Phi ({\varphi }_1^{-1}(\omega ))$. Equivalently, when two normal densities intersect at the threshold value, the probability of misclassifying a genotype 0 subject will not equal the probability of misclassifying a genotype 1 subject, and therefore $li{m}_{n\to \infty }{\hat{p}}_0\ne 0.5=p_0$. For a given threshold, c*, we can define the bias by $p_0^{*}-p_0=p_0{\Phi }_0(c*)+p_1{\Phi }_1(c*)-p_0=p_0({\Phi }_0(c*)-1)+(1-p_0){\Phi }_1(c*)$. Unlike the previous example, with the Mahalanobis distance, Figure 2 shows that describing the bias as a function of ${\sigma }_2^2$ and p₀ can be difficult. In current tests of association, this inequality shows that it possible for $p_0^{n*A}\ne p_0^{n*U}$ even when $p_0^A=p_0^U$, which as we will see, will lead to an inflated type I error and too many low p-values. Start by noting that GWAS test a surrogate hypothesis, $H_0^{*}:p_0^{n*A}=p_0^{n*U}$, not the true hypothesis of interest, $H_0:p_0^A=p_0^U$. Because $p_0^{n*A}$ need not equal $p_0^{n*U}$ when the distributions for cases and controls differ, $H_0^{*}$, the tested hypothesis, can be false even when H₀ is true. Let T* be a standard test statistic for $H_0^{*}$, which is believed to have the following property, $P(T*>t_{\alpha }^{*}|H_0^{*})=\alpha$. Let us make the reasonable assumption that the difference between $P(T*>t_{\alpha }^{*}|H_0,H_0^{*})$ and $P(T*>t_{\alpha }^{*}|H_0^{*})$ is small, or, in words, when $H_0^{*}$ is known to be true, the validity of H₀ has little effect on the distribution of T*. Then,

Here β is a measure of the power to reject $H_0^{*}$ and we assume β > α. Therefore, the current method of rejecting H₀ whenever T* > $t_{\alpha }^{*}$ is actually anti-conservative if the stated p-value is α

Two Genotype Example: Inconsistency

for these scenarios.

Return of Consistency: Modifying the Mahalanobis Distance

With this change, our estimate ${\hat{p}}_0\equiv \frac{1}{n}{\displaystyle {\sum }_{i}1({\hat{G}}_i=0)}$ would now be consistent.

Dependence/Correlation of Called Genotypes

The central limit theorem allows us to be confident that we have an α-level test because $\sqrt{n}({\tilde{p}}_g^q-p_g^q)\approx N(0,p_g^q(1-p_g^q))$ when {G₁,..., G _n is a vector of independent Bernoulli random variables. Again, we have returned to the two genotype scenario to simplify our discussion.

This proof, which uses Jensen's inequality, clearly shows that the two variables, ${\hat{G}}_{i_1}$ and ${\hat{G}}_{i_2}$, are not independent as that would have implied

$P({\hat{G}}_{i_1}={\hat{G}}_{i_2})=P({\hat{G}}_{i_1}=0,{\hat{G}}_{i_2}=0)+P({\hat{G}}_{i_1}=1,{\hat{G}}_{i_2}=1)={(p_0^{n*})}^2+{(p_1^{n*})}^2$. The consequence of this dependence is that $\sqrt{n}({\hat{p}}_g^q-p_g^q)$ does not follow a $N(0,p_g^q(1-p_g^q))$ distribution, and in turn, that $T({\hat{p}}_0^A,{\hat{p}}_0^U)$ neither follows a ${\chi }_2^2$ distribution nor has P($T({\hat{p}}_0^A,{\hat{p}}_0^U)$ > t _α ) = α.

The first term, $nE[var(\frac{{\displaystyle {\sum }_{i=1}^{n}1({\hat{G}}_i=0)}}{n}|C)]$, represents the uncertainty in the true number of subjects with genotype 0, and can be roughly approximated by ${\hat{p}}_0-{({\hat{p}}_0)}^2$. The second term, $nvar(E[\frac{{\displaystyle {\sum }_{i=1}^{n}1({\hat{G}}_i=0)}}{n}|C])$ reflects that there will be a subpopulation, of a random size ~n|Φ ^M (C) - Φ ^M (E[C])|, that is assigned the 'non-ideal' genotype, where a call is labeled 'non-ideal' if it would have been different had the threshold been E [C]. We can approximate this second term, the overall increase in the variance, by nvar(Φ ^M (C)) if P(${\hat{G}}_i$ = 0|C) ≈ Φ ^M (C) and the cor(${\hat{G}}_{i_1}$, ${\hat{G}}_{i_2}$|C) ≈ 0. In experiments, we can estimate nvar(Φ ^M (C)) by bootstrapping samples of C or ${\hat{\Phi }}^M$(C).

Note that ${\hat{p}}_0\equiv {\hat{\Phi }}^M(C)$ and $p_0^{n*}\equiv E[{\hat{\Phi }}^M(C)]$. The appropriate multiple of the first term, n(${\hat{\Phi }}^M$(C) - Φ ^M (C)), should be well approximated by X-E [C], where X ~binomial(n, E [C]). From our own experience, we have seen that the third term, (Φ ^M (E[C]) - E${\hat{\Phi }}^M$(C)]), is a constant close to 0. The second term, (Φ ^M (C) - Φ ^M (E[C])) is the variable which causes deviation from normality. Again, we could approximate the distribution of this term by bootstrapping C.

Next, we offer an example to demonstrate the effect of dependence among called genotypes. Specifically, using the IBD data, we show that a $N(0,E[{\hat{p}}_{AB}](1-E[{\hat{p}}_{AB}]))$ is a poor approximation for the distribution of $\sqrt{n}({\hat{p}}_{AB}-E[{\hat{p}}_{AB}])$. Because we will use real data, we have purposely chosen to discuss $\sqrt{n}({\hat{p}}_{AB}-E[{\hat{p}}_{AB}])/\sqrt{E[{\hat{p}}_{AB}]-{(E[{\hat{p}}_{AB}])}^2}$ instead of T . There are many reasons that T may not follow a ${\chi }_2^2$ distribution, including $p_g^{*}$ being a poor estimate of p _g , the distributions of S _i being far from normal and population substructure. However, for large n, the only reason that $\sqrt{n}({\hat{p}}_{AB}-E[{\hat{p}}_{AB}])/\sqrt{E[{\hat{p}}_{AB}]-{(E[{\hat{p}}_{AB}])}^2}$ will not be approximately normal is if $\{{\hat{G}}_1,\mathrm{...},{\hat{G}}_n\}$ are not independent. Also, we chose to focus on the AB genotype as this is certain to be one of the genotypes with an overlapping cluster.

For each of our 3137 SNPs in the IBD data, we bootstrap 40 samples of 500 subjects, and calculate 40 values of $\sqrt{n}({\hat{p}}_{AB}-E[{\hat{p}}_{AB}])/\sqrt{E[{\hat{p}}_{AB}]-{(E[{\hat{p}}_{AB}])}^2}$, where $\text{E}[{\hat{p}}_{AB}]$ is estimated by $\frac{1}{40}{\displaystyle {\sum }_{k=1}^{40}{\hat{p}}_{AB}}$. The qq-plot in Figure 3 compares these 40 × 3137 values with a N(0,1). The distribution is far from normal, which implies that $\{{\hat{G}}_1,\mathrm{...},{\hat{G}}_n\}$ are dependent. Some SNPs were more likely to contribute skewed values than others, but the top and bottom 100 values (200 total) are from 64 different SNPs, indicating that no one SNP, or small number of SNPs, is responsible for the deviation in the qq plot. In contrast, the qq-plot from well-behaved SNPs, where the contributions to the density of S from the three genotypes were separated, was the expected straight line, showing that it was not the normal approximation skewing the results (Figure 3). Because the magnitude of the observed values were larger than predicted by theory, the practical implication is that tests based on the statistic, T, estimated by $F_{{\chi }_2^2}$(T), will be anti-conservative(i.e. too many significant p-values after adjusting for multiple testing) under the null hypothesis. Here we also note that if S _i were truly normal, the impact of the dependency would be much less. For more details on the origin of dependency, please see Table 1.

Comparing Tests of Association

When genotyping all SNPs and all subjects, the proportion of Bead-Studio 'p-values' below the three thresholds far exceeded 0.005, 0.001, and 0.0005. Even after removing low quality SNPs and allowing missing calls, the proportion of 'p-values' below the three thresholds were 0.015, 0.005, and 0.004. In contrast, LogiCALL eliminated nearly all false positives. The proportion below the three thresholds were 0.004, 0.002, and 0.002. If the majority of these SNPs are presumed to be null associations, logiCALL appears to be the superior method, so long as the power loss is minimal. Assigning conservative 'p-values' to problematic SNPs is nearly equivalent to removing them. However, because of the tendency for there to be multiple, nearly equivalent, maximum likelihood estimates, the relative distances to the AA, AB, and BB genotypes using the single set of maximum likelihood estimates may not be adequate in identifying questionable calls. Therefore, logiCALL gains an advantage by combining two methods for identifying suspect calls. Additionally, it avoids false positives caused by differential bias, where the proportion of missing calls differs between cases and controls. This new method simplifies testing by requiring no preprocessing and testing all SNPs. The power loss from logiCALL depends on the quality of the data. When the statistic for a SNP cleanly separates into an AA, AB, and BB group, there is no power loss. In our IBD example, the 'p-values' reported for rs2066843 and rs2076756, the two SNPs that are believed to be truly associated with IBD were similar for the two methods, Bead-Studio (2.9 × 10^-9 and 5.1 × 10^-10) and logiCALL (1.5 × 10^-8 and 1.6 × 10^-9). Among those subjects meeting the 96% Bead-Studio call rate, logiCALL found no questionable calls.

The likelihood ratio method had mixed results. Clearly, when the distribution of the summary statistic is a mixture of normals, the estimated genotype proportions are asymptotically unbiased. Unfortunately, this method still resulted in an increased number of false positives. However, if we removed those SNPs where at least one call changed when switching from $\hat{\theta }$ to ${\hat{\theta }}_R$, the false positive rate decreased to the expected level. In theory, all calls should be identical and only the resulting likelihoods should differ. When assigning genotypes, the cost, in likelihood, incurred from forcing the vectors of genotype proportions to be equal should be far less than the cost of switching calls. When the reverse is true, and calls switch, the statistic for the three groups cannot be well separated, and the p-value is suspect.