PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of hheKargerHomeAlertsResources
 
Hum Hered. Mar 2012; 73(1): 1–13.
Published online Dec 30, 2011. doi:  10.1159/000332916
PMCID: PMC3268348
A Generalized Sequential Bonferroni Procedure Using Smoothed Weights for Genome-Wide Association Studies Incorporating Information on Hardy-Weinberg Disequilibrium among Cases
Guimin Gao,a* Guolian Kang,b Jiexun Wang,a Wenan Chen,a Huaizen Qin,c Bo Jiang,d Qizhai Li,e Chuanyu Sun,a Nianjun Liu,d Kellie J. Archer,a and David B. Allisond
aDepartment of Biostatistics, Virginia Commonwealth University, Richmond, Va., USA
bDepartment of Epidemiology and Biostatistics, University of Pennsylvania, Philadelphia, Pa., USA
cDepartment of Biostatistics and Epidemiology, Case Western Reserve University, Cleveland, Ohio, USA
dDepartment of Biostatistics, University of Alabama at Birmingham, Birmingham, Ala., USA
eAcademy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, PR China
*Guimin Gao, Department of Biostatistics, Virginia Commonwealth University, PO Box 980032, Richmond, VA 23298-0032 (USA), Tel. +1 804 827 2775, E-Mail ggao3/at/vcu.edu
†Guimin Gao and Guolian Kang contributed equally to this work.
Received June 15, 2011; Accepted September 7, 2011.
Background/Objectives
For genome-wide association studies (GWAS) with case-control designs, one of the most widely used association tests is the Cochran-Armitage (CA) trend test assuming an additive mode of inheritance. The CA trend test often has higher power than other association tests under additive and multiplicative disease models. However, it can have very low power under a recessive disease model in GWAS. Although tests (such as MAX3) robust to different genetic models have been developed, they often have relatively lower power than the CA trend test under additive and multiplicative models. The goal of this study is to propose an efficient method that not only has higher power than the CA trend test under dominant and recessive models but also maintains the power of the CA trend test under additive and multiplicative models.
Methods
We employed the generalized sequential Bonferroni (GSB) procedure of Holm to incorporate information from a Hardy-Weinberg disequilibrium (HWD) test into the CA trend test based on estimating weights from the p values of the HWD test. We proposed to smooth the weights to reduce possible noise.
Results and Conclusions
Results from extensive simulation studies showed that the proposed GSB procedure can achieve the goal described above.
Key Words: Generalized sequential Bonferroni procedure, Genome-wide association studies, Hardy-Weinberg disequilibrium, Multiple testing, Smoothed weights
With the availability of a large number of datasets of genome-wide single-nucleotide polymorphisms (SNPs) for different human diseases, genome-wide association studies (GWAS) have become a useful approach to the identification of genes involved in common human diseases. For GWAS with case-control designs, one of the most widely used association tests is the Cochran-Armitage (CA) trend test [1,2,3], which assumes an additive mode of inheritance. The CA trend test can control the type I error rate when there is departure from Hardy-Weinberg equilibrium (HWE) in the general population from which cases and controls are sampled [3] and can have higher power for additive and multiplicative disease models than other association tests such as the genotype association test using the Pearson's chi-square test statistic. However, this method can have very low power under a recessive disease model (e.g., [4]); the situation becomes even worse in GWAS analysis due to stringent thresholds adjusted for multiple testing.
To detect an association between a marker and a disease under a recessive or dominant disease model, methods based on testing Hardy-Weinberg disequilibrium (HWD) among cases (e.g., [5]) or testing the difference of HWD coefficients between cases and controls [6] have been proposed under the assumption of HWE in the general population. Although these methods have almost no power for testing association under additive or multiplicative disease models, they can provide useful information for detecting association under recessive and dominant disease models. More importantly, this information can be complementary to that provided by the CA trend test. In the development of methods powerful for different disease models, especially for the four commonly used (additive, multiplicative, dominant, and recessive) models, a straightforward idea is to combine the information provided by the CA trend test with that provided by the HWD test into a single test. For example, the Fisher's method [7] has been used to combine p values from the CA trend test and an HWD test [6,8]. Compared to the CA trend test, Fisher's method can have much higher power for a recessive model and relatively higher power for a dominant model. However, it has lower power for additive and multiplicative models. Another example is the weighted average test statistic proposed by Song and Elston [6]. Based on developing an HWD trend test (which compares the HWD coefficient of cases with that of controls), Song and Elston [6] proposed a test statistic which is a weighted average of the CA trend test and the HWD trend test. It seems difficult to estimate the null distribution of the statistic. However, estimating p values by permutations can be computationally prohibitive for GWAS analysis.
Because the true underlying genetic disease models for real data are unknown, methods robust to different models have been developed. A recently widely discussed robust statistic is MAX3 (e.g., [9]). At a biallelic marker, three CA trend test statistics can be calculated based on defining genotypic values assuming recessive, additive, and dominant modes of inheritances have optimal power when the true underlying disease models are recessive, additive, and dominant, respectively [9]. (As stated earlier, the CA trend test statistic assuming additive genotypic values is widely used for GWAS.) The MAX3 of Freidlin et al. [9] is the maximum of the absolute values of the three CA trend test statistics with genotypic values defined under recessive, additive, and dominant models, respectively. Since the asymptotic distribution of MAX3 is very complicated, approximation methods have been developed to estimate p values for MAX3 (e.g., [10,11]). In addition to MAX3, other robust tests, adaptive approaches based on genetic model selection (GMS) and genetic model exclusion (GME), have been proposed under the assumption of HWE in the general population [4,11,12]. These adaptive approaches first select (or exclude) the underlying genetic model by using an HWD trend test of Song and Elston [6] and then apply a CA trend test (or a linear combination of two CA trend tests), which uses genotypic values corresponding to the selected underlying model(s), to analysis of real data. MAX3 and the adaptive approaches have power robust to different genetic models. However, they can be time-consuming when applied to GWAS and can have relatively lower power than the widely used CA trend test statistic (assuming additive genotypic values) when underlying disease models are additive or multiplicative [13], especially for small sample sizes. For example, from the simulation studies of Joo et al. [13] with 500 cases and 500 controls, when the disease locus followed a multiplicative model, the genetic relative risk for the wild homozygote was 1.5, and the minor allele frequency was 0.1, the power of the CA trend test was 0.35, which is 18.6% higher than the power (0.295) of MAX3, in relative pattern.
Because the additive genetic model is often assumed for GWAS in the literature and most identified genetic variants follow an additive mode of inheritance, in the present study, we aimed to develop a computationally efficient method for GWAS that not only has increased power under dominant and recessive disease models in comparison to the CA trend test but also maintains the power of the CA trend test under additive and multiplicative models. Specifically, we estimated a weight for each marker by using the p value of an HWD test using cases only [5], and then we employed the generalized sequential Bonferroni (GSB) procedure of Holm [14] to adjust the p values from the CA trend test by the weights. This step is to achieve higher power than the CA trend test under dominant and recessive models. Further, we proposed to smooth the weights used in the GSB procedure to make the power of the GSB procedure comparable to the CA trend test under additive and multiplicative models. Extensive simulation studies were performed to evaluate the GSB procedure and to compare it with several existing methods. Results from the simulation studies showed that in comparison to the CA trend test, the proposed GSB procedure had much higher power under a recessive model and could have higher power under a dominant model, and that the proposed procedure had power comparable to the CA trend test under additive and multiplicative models. Additionally, the results showed that the GSB procedure was robust to the departure from HWE in the general population. We applied the GSB procedure to the analysis of the genome-wide case-control data for seven complex diseases from the Wellcome Trust Case Control Consortium [15].
2.1. Notation
Herein we use notations similar to those of Zheng and Ng [4]. We consider a genetic case-control design with r cases and s controls genotyped at m independent SNP markers across the whole genome; therefore, n = r + s represents the total number of individuals in the case-control sample. We assume there are no genotyping errors in the sample. For a biallelic SNP marker with alleles A and a, where A is the minor allele, let gi denote a genotype at the marker, where i denotes the number of minor alleles, and (g0, g1, g2) = (aa, Aa, AA). We assume that cases and controls are sampled from a homogeneous general population in which HWE holds and no stratification exists. In this article, f Dgf Cg, and fg denote the frequencies of genotypes g in the case population, control population, and the general population, respectively. Similarly, f DA, f CA, and fA denote the frequencies of allele A in the same three populations, respectively.
2.2 Association Tests for a Single Marker
Here we review two association tests for single marker analysis: the Cochran-Armitage trend test [1,2,3], and the HWD test using cases only [5]. In our proposed method, we used the HWD test among cases instead of the HWD trend test of Song and Elston [6], which compares the HWD coefficients between cases and controls, because our simulation studies suggested that the HWD test using cases only has slightly higher power than the HWD trend test when HWE holds and no stratification exists in the general population.
2.2.1. Cochran-Armitage Trend Test
Zheng and Ng [4] described the CA trend test statistic based on Sasieni [3] as
equation m1
(1)
where (x0, x1, x2) = (0, x, 1) are prespecified values for the three genotypes (aa, Aa, AA); ri and si denote the observed counts of genotype gi in cases and controls, respectively, and ni = ri + si is the number of individuals with genotype gi in the case-control sample. Under the null hypothesis of no association between a marker and the disease, Zx asymptotically follows a standard normal distribution N(0, 1), equation m2and asymptotically follows a χ2 distribution with one degree of freedom. For recessive, additive (or multiplicative), and dominant models, the optimal choices of x are 0, ½, and 1, and the corresponding statistics are denoted as Z0, Z1/2, and Z1, respectively.
Because the underlying disease model is unknown in reality, the additive model assuming additive genotypic values (x = ½) is often used in practice. As stated earlier, in this article, the CA trend test statistic always denotes the statistic assuming additive genotypic values (x = ½) if other types of genotypic values are not arbitrarily assumed. We use equation m3 to denote the CA trend test statistic with x = ½. The statistic Ttrend can have high power when the underlying genetic model is additive or multiplicative, but it has very low power when the underlying genetic model is recessive.
2.2.2. HWD Test Statistic
The HWD test among cases has been used to test the genetic association between markers and a disease (e.g., [5]) when assuming HWE in the general population. For testing HWD at a marker locus among cases, the test statistic [5,16] based on comparing the observed counts of genotypes among cases with the corresponding expected counts can be written as
equation m4
(2)
whereequation m5 is the estimate of the HWD coefficient in cases equation m6 and equation m7, where hat notation denotes the estimate of the variable from the data. Under the null hypothesis of no association between the marker and the disease, if HWE holds in the general population, then ΔD = 0 and THWD asymptotically follows a χ2 distribution with one degree of freedom.
This HWD test can provide valuable information that can be complementary to or different from that provided by the CA trend test for detecting genetic association when the underlying genetic disease model is dominant or recessive [17] but it provides little and no information under the additive and multiplicative disease models, respectively [5]. However, the HWD test statistic THWD in (2) can have an inflated type I error rate when the minor allele frequency is low [18]. To control the type I error rate of the HWD test for a marker with low frequencies, a chi-square statistic corrected for continuity [19] can be used. A better choice is the HWD exact test [16]. Wigginton et al. [18] proposed a fast approach for the HWD exact test that can be applied to GWAS.
2.3. Multiple Testing for Genome-Wide Association Studies
GWAS involve multiple hypothesis testing, i.e., testing hundreds of thousands (or millions) of SNPs simultaneously. Here, we introduce some concepts on multiple hypothesis testing in the context of GWAS. Let p1, p2, …, pm be p values of association tests for m independent markers corresponding to m null hypotheses of H1, H2, …, Hm, where Hi is the null hypothesis of no association between a marker and a disease. The commonly used criterion for error control in multiple testing is the family-wise error rate (FWER), which is the probability of rejecting at least one true null hypothesis (e.g., [20]). The power of a multiple testing procedure can be measured by average power, taken to be the average power of individual tests corresponding to all false null hypotheses [21].
In GWAS, a traditional association test (such as the CA trend test) is often used to test each single marker for its association with a disease. To control the FWER for multiple testing, the significant threshold for each single marker is often adjusted using the Bonferroni procedure. In this study we use the GSB procedure of Holm [14] to control the FWER for GWAS, and more importantly, to adjust p values from the CA trend test by weights calculated from the p values of the HWD test using cases only.
2.3.1. Generalized Sequential Bonferroni Procedure
Here we describe the GSB procedure of Holm [14] in the context of GWAS. Let equation m8 denote the p values for m independent markers calculated by the CA trend test, where the c superscript simply denotes the CA trend test. Given weights w1, w2, …, wm for the tests with hypotheses H1, H2, …, Hm (where we can let wi = 1 or a constant for all i if no weights are provided), define B values as Bi = pci / wi. Let B(1)B(2) ≤ … ≤ B(m) be the ordered B values, with w(1), w(2), …, w(m) being the weights, and H(1), H(2), …, H(m) being the hypotheses corresponding to the ordered B values. The GSB procedure can be described as follows (see also [20]): starting from i = 1, given H(1), H(2), …, H(i –1) have been tested and rejected, if
equation m9
reject H(i); otherwise, accept H(i), H(i + 1), …, H(m), and stop the GSB procedure. The GSB procedure controls the FWER at the nominal significance level of α if the weights are independent of the CA trend tests. From this procedure, we can see that to reject the hypothesis H(i), the corresponding p value from the CA trend test equation m10 must satisfy
equation m11
In other words, if equation m12 irrespective of the size of the corresponding weight w(i), the GSB procedure will not reject the hypothesis H(i).
In the GSB procedure, if the weights are estimated appropriately by use of prior information, the procedure can have higher power than the corresponding Bonferroni procedure [14,20]. Unfortunately, Holm [14] did not provide a method to determine weights for a GSB procedure.
2.3.2. Estimating Weights by Use of Information on HWD Among Cases
Here we propose a method to estimate weights for the GSB procedure in the context of GWAS by using information on HWD among cases. Let equation m13 be the p values from HWD tests for m markers among cases based on equation (2). We define weight wi for the i-th marker as equation m14 and therefore the B value for the i-th marker equation m15. We note that the Fisher's method (online suppl. Appendix A, www.karger.com/doi/10.1159/000332916) is based on considering the null distribution of the logarithm of these B values. We call the GSB procedure with these weights Bi the GSB-HWD procedure.
To control the FWER of the GSB-HWD procedure, we need to prove that the HWD test statistic THWD and the CA trend test statistic Ttrendat a marker are independent under the null hypothesis of no association in order to guarantee that theweights are calculated by using information independent of the CA trend tests. Zheng and Ng [4] proved that when HWE holds in the general population, the two statistics are asymptotically uncorrelated, i.e., the correlation coefficient of Ttrendand THWD equals zero, under the null hypothesis. In this article we further prove that the two statistics Ttrendand THWD are asymptotically independent (see online suppl. Appendix A) by demonstrating the vector (Ttrend, THWD) follows a bivariate normal distribution when HWE holds in the general population.
2.3.3. GSB-HWD Procedure with Smoothed Weights
As will be shown in the simulation study section, in GWAS, the GSB-HWD procedure can have much higher power than the CA trend test (with Bonferroni correction) under the dominant and recessive models. However, like Fisher's method, the GSB-HWD procedure has lower power than CA trend test under the additive and multiplicative models because in this situation the weights from the HWD test provides essentially no information but simply contributes noise to the genetic association studies.
The additive model is often assumed for GWAS in the literature when the true underlying genetic model is unknown. Therefore, we hope to modify the GSB-HWD procedure so that it can maintain the power of the CA trend test under the additive and multiplicative models, but still have increased power under the dominant and recessive disease models. We employ a smoothing method of Roeder et al. [22] to achieve this goal. We calculate a smoothed weight wi for each marker as a linear function of the original weight wi and the average value (w–) of the weights for all markers. Therefore, equation m16 where
equation m17
and λ is a parameter that represents the level of smooth (0 ≤ λ ≤ 1). If λ = 1, then all weights wi equal the constant w–, no weight information from the HWD test is used, and the GSB-HWD procedure is approximately equivalent to the CA trend test and has the highest power for additive and multiplicative models. If λ = 0, then wi = wi, i.e., the weights are not smoothed, and the GSB-HWD procedure will have the highest power for recessive and dominant models. Roughly speaking, the larger the value of λ, the higher the power of the GSB-HWD procedure will have for the additive and multiplicative models and the lower the power for recessive and dominant models. A problem is how to determine the optimal value of λ. Since we hope to preserve the power for additive and multiplicative models, we suggest using λ ≥ 0.6. In this paper, we used λ = 0.6. For simplicity, we call the GSB-HWD procedure using the smoothed weights the S-GSB-HWD procedure.
To evaluate our proposed methods, we compared the FWER and average power of our methods (GSB-HWD procedure and S-GSB-HWD procedure) with those of several existing methods (HWD test, CA trend test, Fisher's method and MAX3) by simulation studies. As in other previous studies [4,23], we used Wright's coefficient of inbreeding F to measure the deviation from HWE in the general population. In our simulation studies, we considered two scenarios: (1) HWE held (F = 0) in the general population, and (2) HWE was slightly violated (we considered F = 0.005 and 0.01).
Given the inbreeding coefficient F and frequency fA of the minor allele A at an SNP marker in the general population, the genotype frequencies fg at the SNP marker in the population can be computed as [23]
equation m18
(3)
We considered two case-control designs: (1) 1,000 cases and 1,000 controls with 1,000 SNPs. This design corresponds to an association study design for candidate genes. (2) 2,000 cases and 3,000 controls with 400,000 SNPs. This design is based on the GWAS data from the Wellcome Trust Case Control Consortium (WTCCC) [15]. Our simulated datasets were generated using a process similar to that described by Joo et al. [13].
3.1. Evaluation of the FWER
To estimate the FWER for the methods being compared, we simulated 100,000 replicated datasets for each case-control design under the assumption that none of the m SNPs was associated with a disease, where m is the number of SNPs. We randomly sampled MAF f for the general population at each of the m SNPs from a uniform distribution U[0.1, 0.50] independently. To generate genotypes at an SNP for r cases and s controls, we followed a process of Joo et al. [13] and randomly generated genotype counts (raa, rAa, rAA) for the cases, where r = raa + rAa + rAA, from the multinomial distribution with parameter r and (faa, fAa, fAA), where faa, fAa, and fAA are calculated by equation (3). We then randomly assign the genotypes aa, Aa, and AA to raa, rAa, and rAA individuals among the r cases, respectively. Similarly, we randomly generated genotype counts (saa, sAa, sAA) and assigned genotypes to the s controls by using the multinomial distribution with parameter s and (faa, fAa, fAA), where s = saa + sAa + sAA. FWERs were estimated as the proportion of 100,000 replicated datasets in which at least one of m true null hypotheses was rejected.
3.2. Evaluation of Power
For each case-control design, we chose 10 SNPs as disease SNPs (i.e., SNPs associated with the disease) with independent genetic effects. We set all disease SNPs to have the same MAFs in the general population. We considered three values for the MAFs: 0.1, 0.3 and 0.5. At each of the remaining non-disease SNPs, we generated genotypes for cases and controls in the same way as described in the previous section (on evaluation of FWER) for generating genotypes under the null hypothesis. Here we describe how we generated genotype data for each disease SNP. Let ψg = Pr(disease [mid ] g) be the penetrance of genotype g (g = AA, Aa, aa) at an SNP marker in the general population. Let γ1 = ψAa / ψaa and γ2 = ψAA / ψaa be the genotype relative risks (GRRs) [24]. For the additive, multiplicative, recessive, and dominant genetic models, γ1 and γ2 satisfy the equations γ2 = 2γ1 − 1, γ2 = γ21, γ1 = 1 and γ1 = γ2, respectively. Therefore, in our simulation studies, we used γ1 as the GRR parameter for the additive, multiplicative and dominant models and γ2 for the recessive model. For each replicated dataset, we set all disease SNPs to have the same GRR values.
Given the population prevalence (k), the inbreeding coefficient (F), population MAF (fA), and GRRs under a disease model, the penetrances ψg can be written as [25]: ψaa = k / (γ2f2A + 2γ1fA (1 – fA) + (1 – fA)2), ψAa = ψaaγ1, ψAA = ψaaγ2. The frequencies of genotype g among cases and controls for a disease SNP can be calculated by equation m19 and equation m20 where fgis genotype frequency in the general population as described above in equation (3). We randomly generated genotypes at each disease SNP marker for r cases and s controls by using the multinomial distributions equation m21 and equation m22, respectively, in a similar way as described in the previous section on evaluating the FWER.
For each case-control design and each set of parameters, we estimated the average power of each method as the average power of individual tests for the 10 disease SNPs based on 1,000 replicated datasets.
3.3. Applications to GWAS of WTCCC Data
The WTCCC [15] conducted GWAS analyses for seven complex human diseases. The dataset for each disease consists of about 2,000 cases and 3,000 shared controls (genotyped with the Affymetrix 500K SNP chip). After quality control filtering, Less than 400,000 SNPs that mapped on the autosomes remained in the dataset for each disease. In this article, we considered only the analysis of autosomal SNPs. By using the threshold α = 5 × 10–7 for each single SNP, the WTCCC identified a total of 21 disease-susceptible regions with significant SNPs for the seven diseases by either the genotype association test (GAT) using Pearson's chi-square test statistic or the CA trend test (see online suppl. table 5). To make comparisons with these results at the same level of threshold, when we applied the GSB procedures to m autosomal SNPs in each dataset, we used a nominal FWER (overall significance) level of 5 × 10–7 · m.
We applied the six methods (HWD test, CA trend test, Fisher's method, GSB-HWD procedure, S-GSB-HWD procedure and MAX3) to the simulated datasets to evaluate FWERs and power. For MAX3, we used the asymptotic distribution approach asy of Zang et al. [11] implemented in the R package Rassoc to estimate the p values.
4.1. Evaluation of FWER
The empirical FWER at the nominal significance levels of 0.01, 0.05, and 0.1, for the two case-control designs based on 100,000 replicates assuming that no SNPs are associated with the disease are reported in table table1.1. When HWE holds (F = 0) or is slightly violated (F = 0.005, 0.01) in the general population, we can see that CA trend test, MAX3 and our two GSB-HWD procedures controlled the FWER quite well. The HWD method had an inflated FWER, especially when HWE did not hold. Fisher's method had a slightly inflated FWER when HWE held and had FWERs much higher than the corresponding nominal levels when HWE did not hold. To control the FWER, our S-GSB-HWD procedure requires only that the weights estimated from the HWD test are independent of the CA trend test, no matter whether the HWD test controls the FWER or not (see also [14]).
Table 1
Table 1
Empirical FWERs based on 100,000 replicates (prevalence k = 0.1)
4.2. Power Comparisons
Since the HWD test using cases only and Fisher's method can have inflated FWERs, especially when HWE does not hold in the general population, we only listed the power of these methods in the tables when HWE holds in the general population to show how much information the HWD test can provide and to show that Fisher's method has lower power than the CA trend test under the additive and multiplicative models.
Results when HWE Held in the General Population
When HWE holds in the general population (F = 0), the estimated average power of the methods being compared is reported in tables tables2,2, ,3,3, ,4,4, ,5.5. From tables tables22 and and3,3, we can see that under the additive and multiplicative models, MAX3 always had lower power than the CA trend test, especially for small sample size, whereas our S-GSB-HWD procedure (using smoothed weights) had power always very close to that of the CA trend test. For example, for the case-control design with 1,000 SNPs, when GRR = 1.25 (i.e., γ1 = 1.25, γ2 = 1.5) and MAF = 0.3, if the disease SNPs followed the additive model, the power of the CA trend test, our S-GSB-HWD, and MAX3 was 0.2732, 0.2612 and 0.2273, respectively (table (table2);2); the power of MAX3 is 16.8% lower than that of CA trend test, in the relative pattern. From table table4,4, we can see that under the recessive model, our S-GSB-HWD procedure had much higher power than the CA trend test; we also noted that MAX3 had higher power than our S-GSB-HWD procedure. For example, when GRR = 1.6 (i.e., γ1 = 1.0, γ2 = 1.6), and MAF = 0.3, for the design with 400,000 SNPs, the power of the CA trend test, our S-GSB-HWD procedure, and MAX3 was 0.1299, 0.4486 and 0.6366, respectively; for the design with 1,000 SNPs, the power of the CA trend test, our S-GSB-HWD and MAX3 was 0.0778, 0.1834 and 0.2817, respectively.
Table 2
Table 2
Empirical average power (for 10 disease SNPs in the data) based on 1,000 replicates (data were simulated under the additive model assuming Hardy-Weinberg equilibrium with prevalence k = 0.1, nominal FWER level α = 0.05)
Table 3
Table 3
Empirical average power (for 10 disease SNPs in the data) based on 1,000 replicates (data were simulated under the multiplicative model assuming Hardy-Weinberg equilibrium with prevalence k = 0.1, nominal FWER level α = 0.05)
Table 4
Table 4
Empirical average power (for 10 disease SNPs in the data) based on 1,000 replicates (data were simulated under the recessive model assuming Hardy-Weinberg equilibrium with prevalence k = 0.1, nominal FWER level α = 0.05)
Table 5
Table 5
Empirical average power (for 10 disease SNPs in the data) based on 1,000 replicates (data were simulated under the dominant model assuming Hardy-Weinberg equilibrium with prevalence k = 0.1, nominal FWER level α = 0.05)
From table table5,5, we can see that when the underlying disease model is dominant, our S-GSB-HWD procedure had power higher than the CA trend test, except when MAF = 0.1, the two methods had comparable power. We can also see that our S-GSB-HWD procedure had power slightly higher than MAX3 (when MAF is equal to either 0.1 or 0.3) and had power lower than MAX3 (when MAF = 0.5 and GRR >1.35). For example, when GRR = 1.3 (i.e., γ1 = γ2 = 1.3), and MAF = 0.3, for the design with 400,000 SNPs, the power of the CA trend test, our S-GSB-HWD procedure and MAX3 was 0.2281, 0.3608 and 0.3384, respectively. As stated earlier, MAX3 had relatively lower power than the CA trend test under the additive and multiplicative models.
Results when HWE Is Violated in the General Population
When HWE was not held in the general population (F = 0.005 and 0.01), the estimated average power of the methods being compared is reported in online supplementary tables 1–4. When F = 0.005 and 0.01, in most situations, the power of our S-GSB-HWD method was very close to the corresponding power of the method with F = 0 (reported in tables tables2,2, ,3,3, ,4,4, ,5),5), except that for the dominant model, the S-GSB-HWD method had decreased power in some cases when F = 0.01 (see online suppl. table 4). However, our S-GSB-HWD method still had power higher than the CA trend test under the recessive and dominant models.
4.3. Application to GWAS of the WTCCC Data
We applied our S-GSB-HWD procedure (see online suppl. table 5) to the analyses of the WTCCC data. Since the HWD test using cases only and Fisher's method cannot control the FWER, and the original GSB-HWD procedure has relatively lower power than the CA trend test under the additive and multiplicative models, we did not apply these three methods to the WTCCC data. The WTCCC group applied the CA trend test and the genotype association test (GAT) using the Pearson's chi-square test statistic (with a threshold 5 × 10–7 for each SNP) to the dataset for each of the seven diseases. The CA trend test identified significant SNPs in 18 regions in total for the seven diseases. All these regions were identified by our S-GSB-HWD procedure using smoothed weights. Our method identified an additional significant SNP rs17042725 for bipolar disease that was not located in the regions identified by WTCCC with the CA trend test or GAT. This SNP may be used for replication or follow-up molecular or bioinformatics studies.
For comparison, we also applied MAX3 to the data (Max was implemented by the asymptotic distribution approach asy of Zang et al. [11]). For Crohn's disease, with the threshold 5 × 10–7 (used by the WTCCC) for single SNPs, the SNP rs10761659 was identified by the CA trend test (p value = 2.677 × 10–14) and our S-GSB-HWD procedure but not by MAX3 (p value = 7.516 × 10–7). This might be because this SNP may follow the additive or multiplicative model, and MAX3 has relatively lower power under the additive and multiplicative models. We also note that some SNPs identified by MAX3 were not identified by the CA trend test and our S-GSB-HWD procedure. In addition, MAX3 took a much longer time than the other methods for the GWAS analysis of the WTCCC data. For example, for the GWAS analysis of the coronary artery disease dataset, MAX3 took about 566 minutes whereas the CA trend test and our S-GSB-HWD procedure each took less than one minute. We used a Beowulf computer cluster [Intel® Xeon® X5550, 2.66 GHz, 8 M cache, turbo, HT, 32 GB memory (8 × 4 GB), 1333 MHz dual ranked RDIMMs for 2 processors, advanced ECC].
We have proposed a simple and fast S-GSB-HWD procedure that uses smoothed weights for GWAS. The S-GSB-HWD procedure is based on application of the generalized sequential Bonferroni procedure of Holm to the incorporation of information of the HWD test using cases only into the CA trend test. When HWE holds or is only slightly violated in the general population, our S-GSB-HWD procedure can control the FWER very well. Compared to the widely used CA trend test (which assumes an additive mode of inheritance), our S-GSB-HWD procedure has much higher power for the recessive disease model and higher or comparable power for the dominant model but also roughly maintains the power of the CA trend test for the additive and multiplicative models. In addition, compared to the CA trend test, for case-control data with larger sample sizes, our S-GSB-HWD procedure can gain much more increased power under the recessive and dominant disease models compared to the CA trend test (tables (tables2,2, ,3,3, ,4,4, ,5).5). Therefore, our proposed method is more suitable for data with large sample sizes (such as ≥1,000 cases and ≥1,000 controls).
In the S-GSB-HWD procedure, the smoothed weights depend on the smooth parameter λ; how to determine the optimal value of λ is an open question. Based on our simulation studies, it appears the optimal value of λ does not depend on the number of SNPs but only depends on the underlying disease models. When the underlying disease models are unknown, setting λ = 0.6 is a good choice.
We note that when the MAF is low (such as ≤10%), in comparison to the CA trend test, our S-GSB-HWD procedure does not have much increased power for all GRR values under the dominant model and for GRR ≤2.0 under the recessive model because the procedure does not gain much information from the HWD test using cases only.
From section 2.3.1.on the GSB procedure, we know that one required condition to reject the null hypothesis of no association at a marker is that the p value from the CA trend test (pc) is less than the nominal FWER level α(such as 0.05), i.e., pc < α (this means that the marker has at least some weak association signal). In other words, if pc ≥ α, no matter how small the p value from the HWD test (or how big the weight for the marker), the S-GSB-HWD cannot reject the null hypothesis. In some situations, genotyping error can cause very small p values of the HWD test using cases only. This feature of the GSB procedure described here can prevent some false positive findings due to the very small p values caused by genotyping errors.
In our S-GSB-HWD procedure, we used p values from the HWD test statistic THWD by using cases only to estimate weights. When HWE holds or is slightly violated, although the HWD test cannot control the FWER for multiple testing, this does not influence that the S-GSB-HWD procedure controls the FWER well. To control the FWER, our S-GSB-HWD procedure requires only that the weights are independent of the p values from the CA trend test. For data with a small sample size or with very small MAF (such as <10%), it might be better to use the HWD exact test statistic [18] to replace the traditional HWD test statistic THWD for more accurate results. When HWE holds in the population, the p values from the exact test may not follow a uniform distribution U[0, 1] under the null hypothesis of no association [18,23] depending on allele frequencies and sample sizes. Therefore these p values cannot be used in Fisher's method, which requires that p values from each test follow U[0, 1] distribution. However, the p values from the HWD exact test can be used to estimate weights in our S-GSB-HWD procedure. We also note that when HWE holds in the general population, if the sample size is small (such as with 1,000 cases), Fisher's method often has slightly violated type I error rate (see our simulation results). When HWE is violated in the general population, Fisher's method has type I error rate much higher than the nominal levels.
In the S-GSB-HWD procedure, we calculate weights on the basis of reciprocal function of p values from the HWD test. We did test other functions, such as –log, square root, and entropy of the p values. It seems that the reciprocal function is optimal because it generated the highest power in our simulation studies. When using the reciprocal function of p values in the S-GSB-HWD procedure, we essentially compare the product of a p value from the HWD test and a p value from the CA trend test with the corresponding threshold. We note that Fisher's method is based on using the distribution of the logarithm of this product. However Fisher's method cannot control the FWER well while our S-GSB-HWD procedure can.
Our S-GSB-HWD procedure is based on the assumptions of independence, i.e., absence of linkage disequilibrium between markers, and cannot handle LD among markers. Wu et al. [26] proposed an SNP-set analysis method which can account for LD and interaction among a set of SNPs and reduce the number of tests in GWAS. Therefore using SNP-set analysis in GWAS can achieve higher power than using individual SNP tests. In our future research, we hope to incorporate the SNP-set analysis into our S-GSB-HWD procedure.
Our S-GSB-HWD procedure also assumes no population stratification and no genotyping error. These issues often exist in the analysis of genetic data for both candidate association studies and GWAS. Since the S-GSB-HWD procedure estimates weights by use of p values from the HWD test using cases only, our S-GSB-HWD procedure may be sensitive to genotyping errors. Song and Elston [6] showed their HWD trend test (detecting the difference of HWD coefficients between cases and controls) can be robust to genotyping errors and population stratification. In addition, Price et al. [27] suggested that a generalization of the CA trend test can be used to control the stratification based on principal components analysis. Our future research will include extending our S-GSB-HWD procedure to handle genotyping errors and stratification. We plan to replace the HWD test using cases only and the CA trend test by the HWD trend test and generalized CA trend test, respectively, in the S-GSB-HWD procedure. However, when there is no obvious stratification in the case-control sample data, especially for data from a homogenous population, use of the HWD test (using cases only) and CA trend test as in our S-GSB-HWD procedure can generate higher power. We note that it may be difficult to extend the (non-permutation-based) approximation methods of Zang et al. [11] for MAX3 to handle population stratification.
Our methods in the present study for GWAS are based on controlling the FWER for multiple testing. It is not difficult to extend the method to control of the false discovery rate.
Here we use the notations defined in Section 2.1. For simplicity, let g [set membership] {0, 1, 2}stands for the genotype that carries g copies of minor allele A at an SNP marker. We define an indicator variable rig = 1 if the i-th case has g copies of allele A and rig = 0 otherwise. Similarly, sg = 1 if the j-th control has g copies of allele A and rg = 0 otherwise. Then the numbers of individuals with genotype g in cases and controls can be calculated as equation m23 and equation m24 respectively. Let equation m25 We assume that HWE holds in the general population. Under the null hypothesis H0 : the marker is not associated with the disease, the hypothesis H01 is true, where H01 is that HWE holds in the cases at the marker.
Apparently, the HWD test stati stic equation m26 (see section 2.2. in the main text), where equation m27 is the estimate of HWD coefficient in the cases equation m28. Then sample variance of DHWD in the cases only equation m29 under H01 : ΔD = 0. Similarly, the CA trend test statistic can be written as equation m30 where equation m31 where equation m32 and equation m33 The estimated variance equation m34 under the null hypothesis H0 because equation m35.
Define equation m36 where (.)'stands for the transpose of the underlying vector and equation m37is the difference between the estimated and actual frequencies of genotype AAin the cases. Then we have
equation m38
(A1)
Under H01 and H0 we observe that
equation m39
(A2)
and
equation m40
(A3)
where
equation m41
(A4)
From the definition, we have equation m42 equation m43 with equation m44 and equation m45 with equation m46 For each g [set membership] {0, 1, 2}, {rig : i = 1, …, r} and {sjg : j = 1, …, s} are two set of independent bounded random variables. Let ξ = (u1, r12, .v1)… Under H01 and H0, we have
equation m47
(A5)
By the multivariate version of the Lindeberg-Levy Central Limit Theorem [28], we obtain the joint asymptotic normality of d:
equation m48
(A6)
Where equation m49 reads as ‘convergences in distribution to’. By the laws of large numbers, we have equation m50 as equation m51 where equation m52 reads as ‘convergences in probability to’. Similarly, we have equation m53 and equation m54
Under H01 and H0, as (r,s)→(∞,∞) we have equation m55 for each genotype g [set membership] {0, 1, 2}, equation m56 Thus, as (r,s)→(∞,∞) and rn−1c under H01 and H0, we obtain that,
equation m57
(A7)
equation m58
(A8)
Substituting (A6)–(A8) into (A3) and applying Slutsky's theorem [28], we obtain the joint asymptotic normality of equation m59 It follows from (A5) and (A7) that ∞ = I2, the 2 × 2 identity matrix.
Supplementary Material
Supplemental Materials
Acknowledgements
This research was supported by grants GM073766, GM077490, and GM081488from the National Institute of General Medical Sciences. Q. Li was partially supported by NSFC No. 10901155. This study makes use of data generated by the Wellcome Trust Case Control Consortium. A full list of the investigators who contributed to the generation of the data is available from www.wtccc.org.uk. Funding for the project was provided by the Wellcome Trust under award 076113 and 085475. We thank Xiangning Chen, Chun-Yu Liu, and Fengjiao Hu for helpful discussion.
1. Cochran WG. Some methods for strengthening the common χ2 tests. Biometrics. 1954;10:417–451.
2. Armitage P. Tests for linear trends in proportions and frequencies. Biometrics. 1955;11:375–386.
3. Sasieni PD. From genotypes to genes: doubling the sample size. Biometrics. 1997;53:1253–1261. [PubMed]
4. Zheng G, Ng HKT. Genetic model selection in two-phase analysis for case-control association studies. Biostatistics. 2008;9:391–399. [PMC free article] [PubMed]
5. Nielsen DM, Ehm MG, Weir BS. Detecting marker-disease association by testing for Hardy-Weinberg disequilibrium at a marker locus. Am J Hum Genet. 1999;63:1531–1540. [PubMed]
6. Song K, Elston RC. A powerful method of combining measures of association and Hardy-Weinberg disequilibrium for fine-mapping in case-control studies. Stat Med. 2006;25:105–126. [PubMed]
7. Fisher RA. Statistical Methods for Research Workers. London: Oliver and Boyd; 1932.
8. Hoh J, Wille A, Ott J. Trimming, weighting, and grouping SNPs in human case-control association studies. Genome Res. 2001;11:2115–2119. [PubMed]
9. Freidlin B, Zheng G, Li Z, Gastwirth JL. Trend tests for case-control studies of genetic markers: power, sample size and robustness. Hum Hered. 2002;53:146–152 (erratum in Hum Hered 2009;68:220). [PubMed]
10. Li Q, Zheng G, Li Z, Yu K. Efficient approximation of P-value of the maximum of correlated tests, with applications to genome-wide association studies. Ann Hum Genet. 2008;72:397–406. [PubMed]
11. Zang Y, Fung WK, Zheng G. Simple algorithms to calculate asymptotic null distributions for robust tests in case-control genetic association studies in R. J Stat Softw. 2010;33 issue 8.
12. Joo J, Kwak M, Zheng G. Improving power for testing genetic association in case-control studies by reducing alternative space. Biometrics. 2009;66:266–276. [PubMed]
13. Joo J, Kwak M, Ahn K, Zheng G. A robust genome-wide scan statistic of the Wellcome Trust Case-Control Consortium. Biometrics. 2009;65:1115–1122. [PubMed]
14. Holm S. A simple sequentially rejective multiple test procedure. Scand J Statist. 1979;6:65–70.
15. The Wellcome Trust Case Control Consortium (WTCCC) Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature. 2007;447:661–683. [PMC free article] [PubMed]
16. Weir BS. Genetic Data Analysis II: Methods for Discrete Population Genetic Data. Sunderland, Mass.: Sinauer Associates, Inc.; 1996.
17. Zang Y, Fung WK, Zheng G. Tail strength to combine two p values: their correlation cannot be ignored. Am J Hum Genet. 2009;84:291–295. [PubMed]
18. Wigginton JE, Cutler DJ, Abecasis GR. A note on exact tests of Hardy-Weinberg equilibrium. Am J Hum Genet. 2005;76:887–893. [PubMed]
19. Yates F. Contingency tables involving small numbers and the χ2 test. J R Stat Soc (suppl) 1934;1:217–235.
20. Kang G, Ye K, Liu N, Allison DB, Gao G. Weighted multiple hypothesis testing procedures. Stat Appl Genet Mol Biol. 2009;8 Article 23. [PMC free article] [PubMed]
21. Roeder K, Wasserman L. Genome-wide significance levels and weighted hypothesis testing. Stat Sci. 2009;24:398–413. [PMC free article] [PubMed]
22. Roeder K, Devlin B, Wasserman L. Improving power in genome-wide association studies: weights tip the scale. Genet Epidemiol. 2007;31:741–747. [PubMed]
23. Rohlfs RV, Weir BS. Distributions of Hardy-Weinberg equilibrium test statistics. Genetics. 2008;180:1609–1616. [PubMed]
24. Schaid DJ, Sommer SS. Genotype relative risks: methods for design and analysis of candidate-gene association studies. Am J Hum Genet. 1993;53:1114–1126. [PubMed]
25. Terwilliger JD, Ott J. Handbook of Human Genetic Linkage. Baltimore: Johns Hopkins University Press; 1994.
26. Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, Lin X. Powerful SNP-set analysis for case-control genome-wide association studies. Am J Hum Genet. 2010;86:929–942. [PubMed]
27. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38:904–909. [PubMed]
28. Serfling RJ. Approximation Theorems of Mathematical Statistics. New York: John Wiley & Sons; 1980.
Articles from Human Heredity are provided here courtesy of
Karger Publishers