Volume 3 Supplement 7
Genetic Analysis Workshop 16
Genomewide association studies of rheumatoid arthritis data via multiple hypothesis testing methods for correlated tests
 Guolian Kang^{1},
 Douglas K Childers^{1},
 Nianjun Liu^{1},
 Kui Zhang^{1} and
 Guimin Gao^{1}Email author
DOI: 10.1186/175365613S7S38
© Kang et al; licensee BioMed Central Ltd. 2009
Published: 15 December 2009
Abstract
Genomewide association studies often involve testing hundreds of thousands of singlenucleotide polymorphisms (SNPs). These tests may be highly correlated because of linkage disequilibrium among SNPs. Multiple testing correction ignoring the correlation among markers, as is done in the Bonferroni procedure, can cause loss of power. Several multiple testing adjustment methods accounting for correlations among tests have been developed and have shown improved power compared to the Bonferroni procedure. These methods include a Monte Carlo (MC) method and a method of computing pvalues adjusted for correlated tests. The objective of this study is to apply these two multiple testing methods to genomewide association study of the Genetic Analysis Workshop 16 rheumatoid arthritis data from the North American Rheumatoid Arthritis Consortium, to compare the performance of these two methods to the Bonferroni procedure in identifying susceptibility loci underlying rheumatoid arthritis, and to discuss the strengths and weaknesses of these methods. The results show that both the MC method and pvalues adjusted for correlated tests method identified more significant SNPs, thus potentially have higher power than the corresponding Bonferroni methods using the same test statistics as in the MC method and pvalues adjusted for correlated tests, respectively. Simulation studies demonstrate that the MC method may have slightly higher power than the pvalues adjusted for correlated tests method.
Background
Genomewide association studies (GWAS) for complex diseases involve multiple hypothesis testing. The Bonferroni procedure is commonly used to control familywise error rate (FWER) for multiple hypothesis testing. However, the Bonferroni procedure becomes more conservative as the number of hypotheses tested increases and the test statistics are correlated [1, 2]. To handle the correlation among test statistics, a permutation method [3] was proposed based on estimation of the joint distribution of test statistics. However, this approach is computationally intensive and not appropriate for GWAS. Therefore, several efficient methods that can account for the correlation among test statistics have been developed for multiple testing [1, 2]. Lin [2] proposed a Monte Carlo (MC) sampling approach based on approximating the joint distribution of test statistics. This method does not require repeated analyses of simulated datasets as in the permutation method, and therefore is much less computationally demanding. Conneely and Boehnke [1] proposed a method of computing pvalues adjusted for correlated tests (p_ACT) by numerical integration of the asymptotic multivariate normal distribution of the test statistics. This approach is very computationally efficient and attains even greater speed. In this study we applied three multiple testing procedures, the Bonferroni procedure, MC method, and p_ACT method, to GWAS of the Genetic Analysis Workshop 16 (GAW16) rheumatoid arthritis (RA) data from the North American Rheumatoid Arthritis Consortium (NARAC). We compared the performance of these three procedures by simulation studies.
Methods
We describe the three multiple hypothesis testing procedures in the context of association studies. Suppose there are n individuals with m markers in the observed casecontrol data. We test m null hypotheses H_{1}, H_{2}, ..., H_{ m }for the m markers. The corresponding pvalues are p_{1}, p_{2}, ..., p_{ m }. In the Bonferroni procedure, if p_{ i }≤ α/m, then H_{ i }is rejected (i = 1, ..., m), where α is the preset significance level. While in the Bonferroni procedure all tests are assumed to be independent, the MC method and the p_ACT methods described below account for dependence among test statistics by considering the joint distribution of test statistics. All the three methods can control the FWER well.
MC method
where , , Y_{ i }is the phenotypic value of individual i and ; X_{ ji }is the genotypic score of individual i at locus j and ; . If the hypothesis H_{ j }is true, the statistic T_{ j }has approximately a χ^{2} distribution with d_{ j }degrees of freedom, where d_{ j }is the dimension of U_{ j }. For GWAS of the RA data in this article, we only consider an additive genetic model with d_{ j }= 1, and X_{ ji }= 0, 1, or 2, indicating the number of minor alleles in the genotype of individual i at locus j.
The test statistics (T_{1}, T_{2}, ..., T_{ m }) may be correlated due to linkage disequilibrium (LD) among markers. The multiple testing procedure using the actual joint distribution of (T_{1}, T_{2}, ..., T_{ m }) can be computationally intensive. The MC method provides an approach to approximate the actual joint distribution by MC sampling. The MC method defines , where and G_{1}, G_{2}, ..., G_{ n }are independent standard normal random variables that are independent of the data, and then the method uses the joint distribution of values to approximate the joint distribution of T_{ j }values based on obtaining realizations from distributions of values by repeatedly generating the normal random samples G_{1}, G_{2}, ..., G_{ n }. Let t_{(1)} ≥ t_{(2)} ≥ ⋯ ≥ t_{(m) }be the ordered observed values of the test statistics (T_{1}, T_{2}, ..., T_{ m }), and let H_{(1)}, H_{(2)}..., H_{(m) }be the hypotheses and are variables corresponding to (t_{(1)}, t_{(2)}, ..., t_{(m)}), respectively. The MC method works as a stepdown procedure as follows: starting with hypothesis H_{(1)}, the method rejects H_{(j)}, j = 1, 2, ..., m and removes the corresponding marker and variable from consideration, if , provided that H_{(1)}, ..., H_{(j1) }have been tested and rejected. This probability is calculated based on a large number (e.g., 10,000) of realizations of the values. We have implemented this method by using the statistical package R [4].
p_ACT method
where , . These test statistics (T_{1}, T_{1}, ..., T_{ m }) asymptotically follow multivariate normal distribution N(0, R), where R = {r_{ jk }}, j = 1, ..., m, k = 1, ..., m, and ; . The p_ACT method also works as a stepdown procedure: 1) If p_{ ACT }<α, then reject the null hypothesis associated with p_{min}, and remove p_{min} and the marker associated with p_{min}. 2) Let p_{min} = p_{(2)} in Eq. (2), change m into m1, and repeat Step 1. Continuing in this fashion for the remaining p_{(j) }until for some p_{(k)}, p_{ ACT }= α, then accept all remaining hypotheses (including that associated with p_{(k)}). This p_ACT method has been implemented in a computer program in R [6].
Partitioning of genomewide singlenucleotide polymorphism (SNP) data
As stated earlier, the p_ACT method can only handle up to 1,000 tests at a time, and the MC method can handle more than 1,000 tests but may become computationally intensive when the number of tests is very large. To apply these methods to GWAS of the RA data, we divided the whole genome into small blocks; each block includes hundreds or up to one thousand of SNPs. We assume that tests within each block are dependent, and that tests from different blocks are independent. To control the FWER at α for the whole genome, we apply Bonferroni procedure among blocks, that is, we assign α_{ b }to each block such that the sum of these α_{ b }equals α (i.e., ∑ α_{ b }= α), where α_{ b }is proportional to the number of SNPs within the block (i.e., block size). We applied the MC method and p_ACT method separately to each block and control FWER at α_{ b }for the tests within the block. In this study, we considered the block sizes of 100, 500, and 1,000 separately to evaluate the effect of block size on the performance of the MC method and p_ACT method.
Application to RA data
The RA dataset contains 868 cases and 1,194 controls with 545,080 SNPs after removing duplicated and contaminated samples. If an individual has missing genotype at a marker, we imputed the mostfrequent genotype observed in the data at that marker. We removed SNPs with minor allele frequencies (MAF) less than 0.01 and SNPs with pvalues less than 1 × 10^{4} in HardyWeinberg equilibrium (HWE) test in controls. We did not consider SNPs on sex chromosomes. After these procedures, 515,050 SNPs remained in our analysis.
We applied the three multiple testing procedures to the GWAS of the RA dataset. First, we performed the association analysis on the RA dataset by using the test statistics in Eqs. (1) and (3) separately and obtained pvalues for each SNP. We call these pvalues raw pvalues. For each test statistic, we applied the Bonferroni procedure to the pvalues and set the nominal level of FWER as 0.05. The raw pvalues calculated from statistic in Eq. (3) were used in the p_ACT method (see below). The MC method is based on the statistic in Eq. (1), and we set the number of replicates of the normal random samples G_{1}, G_{2}, ..., G_{ n }as 25,000. The p_ACT method is based on the test statistic in Eq. (3), and all tests are twosided. In the p_ACT method we set the limit on the number of simulations or integrand values in R function "pmvnorm" as maxpts = 25,000.
Simulation studies
To evaluate the performance of the two statistics in Eqs. (1) and (3) and of three multiple testing methods, we simulated 10,000 replicated data sets in a manner similar to Lin [2]; each data set included N_{1} cases and N_{2} controls with 100 SNPs in a chromosomal region (or one block). We considered two sets of values of N_{1} and N_{2}: 1) N_{1} = N_{2} = 100, and 2) N_{1} = 100, and N_{2} = 150. For each individual, the simulated chromosomal region consisted five independent consecutive subregions. Each subregion has 20 biallelic SNPs in LD with coefficient of r^{2} = 0.9 between two successive SNPs. At each SNP, HWE was assumed and MAF was 0.3. We chose one SNP in the first subregion as a disease SNP, and determined case/control status based on an additive disease model with disease prevalence of 0.1 and genotype relative risk of 2.
Results
Estimated FWER and power from the simulated 10,000 replicated data sets
The numbers of identified significant SNPs from the RA data set
We compared computing times of the MC method and p_ACT method. As an example, we only showed the times for chromosome 9. With block sizes of 100, 500, and 1,000, the MC method used about 6.24 hr, 2.89 hr, and 2.48 hr, while the p_ACT method used 0.27 hr, 0.47 hr, and 1.08 hr, respectively. The p_ACT method is faster than the MC method.
Discussion
We have applied two multiple testing methods (MC and p_ACT), which account for correlation among tests by splitting each chromosome into smaller blocks, to the GWAS of the RA dataset and then compared the results of these methods to those of the Bonferroni procedure. Both the MC method and p_ACT method identified more significant SNPs than the Bonferroni procedure. The numbers of significant SNPs identified by the MC and p_ACT methods increased as the block size increased.
The test statistic in Eq. (3) is transformed from a traditional score statistic from a generalized linear model. The essential difference between the statistics in Eqs. (1) and (3) is that variance V_{ j }is estimated by different ways. Our simulation studies show that when the numbers of cases and controls are equal, the two test statistics have almost the same power, and that when the numbers of cases and controls are different, the test statistic in Eq. (1) can have slightly higher power than that in Eq. (3), and consequently, the power of the MC method can be slightly higher than that of p_ACT method. Our simulation studies were only based on additive model, small sample sizes, and small number of SNPs. More extensive simulation studies are necessary in the future research.
In our analysis, we divided the whole genome into blocks with a fixed number of SNPs. We only accounted for LD within each block, and we assumed independence between tests from different blocks by ignoring LD between blocks. This assumption can cause loss of power. To avoid losing power, each entire chromosome may be treated as a block. However, the MC method will become computationally intensive or infeasible, and the p_ACT cannot handle more than 1,000 SNPs in each block. Another possible solution is to split each chromosome into blocks according to LD pattern, to group a set of consecutive SNPs in strong LD into one block and to ignore weak LD between blocks. As described earlier, the larger the block size we select, the higher the power we can obtain. This is an issue we will pursue in the future. Also we did not consider population stratification, which may cause spurious falsepositive results.
List of abbreviations used
 FWER:

Familywise error rate
 GAW16:

Genetic Analysis Workshop 16
 GWAS:

genomewide association studies
 HWE:

HardyWeinberg equilibrium
 LD:

Linkage disequilibrium
 MAF:

Minor allele frequency
 MC:

Monte Carlo
 NARAC:

North American Rheumatoid Arthritis Consortium
 p_ACT p:

Values adjusted for correlated tests
 RA:

Rheumatoid arthritis
 SNP:

Singlenucleotide polymorphism.
Declarations
Acknowledgements
This research was supported by grants GM073766, GM074913, GM081488, and GM077490 from the National Institute of General Medical Sciences and T32HL072757 from the National Heart, Lung, and Blood Institute.
This article has been published as part of BMC Proceedings Volume 3 Supplement 7, 2009: Genetic Analysis Workshop 16. The full contents of the supplement are available online at http://www.biomedcentral.com/17536561/3?issue=S7.
Authors’ Affiliations
References
 Conneely KN, Boehnke M: So many correlated tests, so little time! Rapid adjustment of pvalues for multiple correlated tests. Am J Hum Genet. 2007, 81: 11581168. 10.1086/522036.PubMed CentralView ArticlePubMedGoogle Scholar
 Lin DY: An efficient Monte Carlo approach to assessing statistical significance in genomic studies. Bioinformatics. 2005, 21: 781787. 10.1093/bioinformatics/bti053.View ArticlePubMedGoogle Scholar
 Westfall PH, Young SS: Resamplingbased Multiple Testing: Examples and Methods for pValue Adjustment. 1993, New York, WileyGoogle Scholar
 The Comprehensive R Archive Network. [http://cran.rproject.org]
 Genz A, Bretz F, Hothorn T: mvtnorm: multivariate normal and t distribution. R package version 0.80. [http://cran.rproject.org/doc/packages/mvtnorm/mvtnorm.pdf]
 p_ACT. [http://csg.sph.umich.edu/boehnke/p_act.php]
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.