Power analysis for genome-wide association studies

Abstract

Background

Genome-wide association studies are a promising new tool for deciphering the genetics of complex diseases. To choose the proper sample size and genotyping platform for such studies, power calculations that take into account genetic model, tag SNP selection, and the population of interest are required.

Results

The power of genome-wide association studies can be computed using a set of tag SNPs and a large number of genotyped SNPs in a representative population, such as available through the HapMap project. As expected, power increases with increasing sample size and effect size. Power also depends on the tag SNPs selected. In some cases, more power is obtained by genotyping more individuals at fewer SNPs than fewer individuals at more SNPs.

Conclusion

Genome-wide association studies should be designed thoughtfully, with the choice of genotyping platform and sample size being determined from careful power calculations.

Background

One goal of modern human genetics is to identify the genetic variants that predispose individuals to develop common, complex diseases. It has been proposed that population-based association studies will be more powerful than traditional family-based linkage methods in identifying such high-frequency, low-penetrance alleles [1]. Such studies require the genotypes a large number of polymorphisms (usually single nucleotide polymorphisms [SNPs]) across the genome, each of which is tested for association with the phenotype of interest. As originally proposed, this would be a direct test of association, in which the functional mutation is presumed to be genotyped. An alternate approach to association studies takes advantage of the correlation between SNPs, called linkage disequilibrium (LD), that can occur due to the genealogical history of the polymorphisms [2]. In this approach, often called indirect association, one SNP is genotyped and used to infer indirectly the genotypes at other SNPs with which it is in high LD [3]. As one genotyped SNP, called a "tag" SNP, can be in LD with numerous other SNPs, much fewer SNPs (10⁵ – 10⁶) would need to be genotyped to capture the common variation in the genome [3]. Recent advances in genotyping technology make such studies feasible [4, 5] and the first results of such studies are being published [6–10].

One key question in designing such studies is the choice of tag SNPs. Numerous methods for choosing the best set of tagging SNPs have been developed and compared [11]. One common measure evaluates the pairwise LD, measured by r², between the tag SNPs and all other SNPs [12]. The value r² represents the correlation between two SNPs. It is a useful measure because, if N individuals are needed for a specific power with a direct test of association, N/r² individuals would be needed for an indirect test of association [2]. Sets of tag SNPs are usually compared by their "coverage," or fraction of variants in the genome that are in LD (r² above some threshold) with at least one tag [12–15].

There are two related problems with this measure of coverage. First, the binary decision of whether r² is above or below a threshold does not capture the continual decrease in power as r² decreases. If the cutoff value of r² is 0.8, a SNP that shows LD of r² = 0.75 with a tag would be called undetectable since the measure of LD is below the threshold. In truth, association would be detectable, albeit with reduced power. Second, knowledge of the coverage of a set of tag SNPs says nothing about the number of individuals needed for a well-powered study. A better measure to evaluate tag SNPs would be an explicit calculation of the probability that a genome-wide association study will find a statistically significant association given that such an association exists (i.e., power). To solve this problem, one needs to be able to calculate the power of a study given a specified genetic model and sample size. Skol et al. have proposed a method for computing power, though they were concerned with issues of study design rather than tag SNP choice [16]. Jorgenson and Witte, who noted the same problems, propose a "cumulative r² adjusted power" that integrates LD and tag SNP information to provide the overall power of a study [17].

Realistically, one does not have an unlimited choice of SNPs but rather chooses among several competing commercial products with fixed sets of tag SNPs. Therefore, instead of choosing a set of tag SNPs, a more common problem now is how to evaluate which of several fixed sets of tag SNPs is better for a particular study. Several papers have looked at power for hypothetical and commercial sets of tag SNPs through empirical simulations on a subset of chromosomal regions [13, 18]. This approach suffers from both the speed problem of empirical simulations and the assumption that the sampled regions are representative of the genome as a whole. What is needed is an application of explicit power calculation methods (such as that of Jorgenson and Witte [17]) to the commercially available sets of tag SNPs to allow comparison among products and power calculation for real studies.

Here, I present a method for computing the power of a genome-wide association study when a genetic model and sample size are specified and LD information is available for the population being studied. This method is equivalent to the cumulative r² adjusted power of Jorgenson and Witte [17], which will be referred to as "power" for brevity. I show that to obtain the best power, different commercial genotyping products should be used for different populations. I further find that power is sometimes improved by genotyping more individuals at fewer SNPs rather than fewer individuals at more SNPs. These calculations can guide the optimal design of future genome-wide association studies.

Results and discussion

The power calculations require genotype data on a large representative sample of common SNPs from the population as well as a list of which of these representative SNPs are the tag SNPs (SNPs to be genotyped). Power is computed in three steps. First the best tag SNP for each of the representative SNPs is found. Then, the power for detecting association for each of the representative SNPs assuming that SNP directly influences the phenotype is computed. For this computation, it is assumed that the study will be performed by testing for genotype frequency differences between cases and controls using a two-degree of freedom χ² test in which multiple tests are corrected for using the Bonferroni correction. This test explicitly assumes a codominant model. I use this test because it is the most general, at the cost of reduced power relative to a model-specific test. While a multimarker tagging approach could be taken [13], this added level of complexity is not usually included in a first-pass analysis of genome-wide association data and is therefore including it in our power-calculation would inflate the power one might expect in real-world application of genome-wide association studies. Finally, the average power over all the SNPs is taken to be the power of the study.

Taking the average power over all the SNPs is justified using probability theory. Assume there are N SNPs present in a given population, each one represented as S _i . Let C _i represent SNP i being causative, and D _i represent SNP i being detected. Assume that one of these SNPs is the causative SNP, but it is unknown which of these is the causative SNP. Then the overall power of the study is given by {\displaystyle \sum _{i=1}^N\mathrm{Pr}(C_i,D_i)}. The power computed for a specific SNP S _i is given as P _i = Pr(D _i |C _i ). Thus, if each P _i multiplied by Pr(C _i ), we get

\begin{array}{l}\mathrm{Pr}(C_i,D_i)=\mathrm{Pr}(D_i|C_i)\mathrm{Pr}(C_i)=P_i\mathrm{Pr}(C_i)\hfill \\ Power={\displaystyle \sum _{i=1}^N{P}_i\mathrm{Pr}(C_i)}\hfill \end{array}

The added assumption that each SNP is equally likely to be causative yields

\begin{array}{l}\mathrm{Pr}(C_i)=\frac{1}{N}\hfill \\ Power={\displaystyle \sum _{i=1}^N{P}_i\frac{1}{N}}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{P}_i}\hfill \end{array}

This final equation is the same as taking the average power over all the SNPs.

This method was applied to examine the power of genome-wide association studies in the four populations studied in the International HapMap Project [19]. I examined the performance of the tag SNPs provided by the major high-density genotyping platforms available commercially: 100 K and 500 K SNP sets from Affymetrix and 300 K and 550 K SNP sets from Illumina. (Since then, more products have come on the market; the same approach can be taken with them.) I first asked how many SNPs on each of these arrays would be useful for studying a given population by asking what percentage of tag SNPs provided by each platform are common (minor allele frequency > 5%) in each of the four HapMap populations (Table 1). The largest fraction of common SNPs is found when the Illumina chip is used in the CEU population. As the Illumina chip was designed to optimize coverage of the CEU population, this result is unsurprising.

Table 1 The number of SNPs present in each population and present in each commercial genotyping system

I next asked how power changes with increasing sample size for the various genotyping platforms (Figure 1), populations, and models. For all sets of tag SNPs, as expected, power increases both as the sample size increases and as the magnitude of effect, measured by the genotype relative risk (GRR), increases. While Figure 1 only shows this data for a multiplicative model in the CEU population, similarly shaped curves were observed in the other populations and for other models [see Additional file 1]. In the Affymetrix 500 K and Illumina 300 K SNP sets, the slope of the power curve starts leveling off (approaching zero) with a few thousand individuals when GRR is more than 1.5. For smaller GRRs, the sample sizes required for adequate (at least 50%) power becomes quite large.

Figure 1

Power for the test of genotypic association as a function of sample size at different genotype relative risks (GRR). All panels are for the CEU HapMap population when the number of cases equals the number of controls and a multiplicative model is used. (A) Power for the Affymetrix 100 K system. (B) Power for the Illumina 300 K system. (C) Power for the Affymetrix 500 K system. (D) Power for the Illumina 550 K system.

One critique of this approach is that the non-specific test used may not be the most powerful approach if we know the genetic model the disease follows. For instance, to study a trait that we believe follows a multiplicative model; a 2 × 2 contingency table to test for allelic association may be more appropriate. Power calculations for this test (Figure 2) shows that the relative pattern is the same as for a test of genotypic association, but the power is generally increased when an allelic test is used in instead of a genotypic test. Similar power calculations can be done if one wants to use an explicit test for a dominant or recessive mode of inheritance. However, as can be seen in this comparison between the Affymetrix 500 K and Illumina 550 K genotyping system, choice of SNPs and sample size can play a bigger role in determining power than choice of test. For the specified GRR of 1.5, the Illumina 550 K system with a genotypic test is more powerful than the Affymetrix 500 K system when sample size is greater than 2000 individuals (Figure 2).

Figure 2

Power for genotypic and allelic tests. Data is shown for a GRR of 1.5 under a multiplicative model, the CEU HapMap population, and the specified genotyping system.

Another possible criticism of this method is that the SNPs genotyped as part of the International HapMap Project may not be a representative subset of the common SNPs in the genome as a whole. To investigate this possibility, I compared the coverage of the various SNPs in the ENCODE and non-ENCODE regions from the HapMap project (Figure 3). Since the ENCODE regions of the HapMap project were completely resequenced in a subset of 48 individuals, I hypothesized that almost all common (minor allele frequency >5%) variants would have been identified in that region. If the SNPs genotyped as part of the HapMap are a representative subset of all of the common SNPs, then the coverage of an arbitrary set of tag SNPs should be equal for the two data sets. Assuming tag SNPs were chosen similarly for the ENCODE and non-ENCODE regions, relying on the HapMap data slightly overestimates r² with the tag SNPs and therefore could slightly inflate the power estimation. As the fraction of SNPs with an r² greater than the cutoff differs between the ENCODE and non-ENCODE regions by at most ten percentage points, and an average of three percentage points, this overestimation is not likely to be extreme.

Figure 3

Coverage of tag SNPs. Fraction of non-tag SNPs in LD with a tag SNP with r² above specified threshold for the ENCODE and non-ENCODE regions of the HapMap project for the CEU and YRI populations. Results are shown for the Illumina 550 K (A) and Affymetrix 500 K (B) chips. The JPT+CHB population was not included because the curves generally overlap with the CEU curves and would make the graph harder to read. Results for the JPT+CHB population and for the other chips are qualitatively similar to the curves shown here.

An easy and useful way to compare the power of different tag SNP sets in different populations is the sample size needed to achieve 80% power. The Illumina 550 K clearly performs best in all three populations (Figure 4). For the CEU population, the Illumina 300 K outperforms the Affymetrix 500 K, while in the other two populations the Affymetrix 500 K is better. This is not surprising, as the Illumina chips were optimized on CEU HapMap data. As the Affymetrix 500 K set is really two independent 250 K sets, I also looked at the power of each 250 K set individually. While the complete 500 K set of SNPs has more power than either half, the number of individuals required for 80% power using one half of the set is never twice the number required for the full set. This means that in cases when the number of chips that can be run rather than number of available samples is the limiting factor, it might make more sense to genotype more individuals using only one chip than to genotype fewer individuals using both chips. To test this hypothesis, I plotted power versus the number chips needed for the components of the Affymetrix 500 K system (Figure 5). The number of chips is simply the sample size for Nsp and Sty alone, and twice the sample size for the Nsp+Sty combination. Except in cases where power gets very high due to a large GRR and/or sample size, for a constant number of chips using only one of Nsp or Sty on more individuals provides a more powerful study.

Figure 4

Total individuals required for 80% power. The computations assume the number of cases equals the number of controls and a GRR of 1.75. CEU, JPT+CHB, and YRI are the HapMap populations. Affy 250 K Nsp and Affy 250 K Sty represent the two chips that make up the Affymetrix 500 K genotyping system.

Figure 5

Power as a function of number of chips needed for the Affymetrix 500 K system and its two components. Calculations are done for a GRR of (A) 1.5 and (B) 2.0.

I have presented a method to compute the power of a genome-wide association study in which a fixed set of tag SNPs will be genotyped. For the sake of simplicity, I only considered one straightforward single-SNP analysis scheme. While this approach has been used successfully [6], others have suggested that greater power can be obtained by looking at multiple tags or haplotypes [18, 20]. This method for computing power can be adapted to such strategies provided it is possible to compute the power of detecting each SNP in the population given the set of tagging SNPs. I also assume that each SNP is equally likely to be functional. If we knew a priori the probability that a given SNP is functional, we could use this to weight the average power over all the SNPs. Such a weighting scheme would prioritize SNPs more likely to be of interest because of either functional considerations or location [21]. For instance, assume we assigned each SNP a probability of being the causative SNP based on external evidence such as a prior linkage study. If these probabilities are normalized to sum to one, they can be used to compute a weighted average power in this approach.

Conclusion

Proper design of a genome-wide association study requires careful calculation of the power. These calculations will be invaluable to anyone who is planning a genome-wide association study. Using these calculations, the proper sample size to get adequate power in a given study can be computed. Furthermore, the performance of different genotyping platforms can be compared, allowing an investigator to choose whatever is best for his or her study. By performing such calculations, genome-wide association studies can be optimized to get the maximal power possible for a given set of resources.

Methods

Genotype data and populations

I used genotype data from release 21 (phase II) of the International HapMap project [19]. I used data from all four populations studied in the HapMap project. These populations are defined by the HapMap project as follows: Yoruba in Ibadan, Nigeria (abbreviation: YRI); Japanese in Tokyo, Japan (abbreviation: JPT); Han Chinese in Beijing, China (abbreviation: CHB); and CEPH (Utah residents with ancestry from northern and western Europe) (abbreviation: CEU). Similar to the analysis performed by the HapMap project, I combined genotypes from the JPT and CHB populations to make a joint JPT+CHB population. For all three resulting populations, I removed SNPs that have a minor allele frequency (MAF) less than 0.05 in that population. The remaining SNPs are considered to be "common." A summary of the number of SNPs remaining for each population is found in Table 1. When phased data is needed, I used the phased chromosomes for release 21.

Calculation of power

To compute the overall power of an association study, I use three steps. First, I find the best tag SNP for each genotyped SNP in the data set. Then, I compute the power for each SNP assuming the specified GRR and sample size. Finally, I take an average power over all the SNPs to get the overall power.

To find the best tag SNP for each genotyped SNP, I look at the linkage disequilibrium between each SNP and all tag SNPs within 300 kb of it. For each pair of SNPs, I infer the two-locus haplotype frequencies between them using expectation maximization and compute r² between the two SNPs from the inferred haplotype frequencies [12]. The best tag is then taken to be the tag SNP with the highest value of r².

To compute the power for a SNP, I assume that we are looking at genotype frequency differences using a two-degree of freedom χ² test. The power of this test is computed using a non-central χ² distribution with non-centrality parameter λ. Equations for λ have been previously derived for a general χ² test [22] and for application to genetic association [23]. Specifically, for genotypic association λ is given by:

\lambda =N_A{N}_U[\frac{{(p_{00}-p_{10})}^2}{N_A{p}_{00}+N_U{p}_{10}}+\frac{{(p_{01}-p_{11})}^2}{N_A{p}_{01}+N_U{p}_{11}}+\frac{{(p_{02}-p_{12})}^2}{N_A{p}_{02}+N_U{p}_{12}}]

where N _A and N _U are the number of case (affected) and control (unaffected) individuals, respectively; p₀₀, p₀₁, and p₀₂ are the genotype frequencies in the cases; and p₁₀, p₁₁, and p₁₂ are the genotype frequencies in the controls. If, instead of a 3 × 2 table we use a 2 × 2 table for a one-degree of freedom test of allelic association, the non-centrality parameter is given by

\lambda =2N_A{N}_U{(p_A-p_U)}^2\frac{N_A+N_U}{(N_A{p}_A+N_U{p}_U)(N_A+N_U-N_A{p}_A-N_U{p}_U)}

where p _A and p _U are the frequencies of allele 0 in the cases and controls, respectively.

I use the Bonferroni correction for multiple testing and require a p-value of 0.05/M (where M is the number of tag SNPs genotyped) for statistical significance. When association is directly tested (the SNP is a tag SNP), I use the actual number of cases and controls to compute the power. For indirect association (the SNP is in LD with a tag SNP), I reduce the number of cases and controls by a factor of r² for the power computation [2].

I assume that the disease has a low enough prevalence in the population that the risk allele frequency in those without the disease approximates the risk allele frequency in the population. I can set the disease to follow a multiplicative, additive, dominant, or recessive model with a specified genotype relative risk (GRR) for the SNP of interest [1]. Given that genotype 0 is the wildtype, and taking p₁₀, p₁₁, and p₁₂ from the observed genotype frequencies in the population, p₀₀, p₀₁, and p₀₂ are computed as follows:

Multiplicative

\begin{array}{l}{p}_{00}=\frac{\frac{p_{10}^2}{p_{11}{p}_{12}}}{{\gamma }^2\frac{p_{10}}{p_{11}}+\gamma \frac{p_{10}}{p_{12}}+\frac{p_{10}^2}{p_{11}{p}_{12}}}\hfill \\ p_{01}=\frac{\gamma \frac{p_{10}}{p_{12}}}{{\gamma }^2\frac{p_{10}}{p_{11}}+\gamma \frac{p_{10}}{p_{12}}+\frac{p_{10}^2}{p_{11}{p}_{12}}}\hfill \\ p_{02}=\frac{{\gamma }^2\frac{p_{10}}{p_{11}}}{{\gamma }^2\frac{p_{10}}{p_{11}}+\gamma \frac{p_{10}}{p_{12}}+\frac{p_{10}^2}{p_{11}{p}_{12}}}\hfill \end{array}

Additive

\begin{array}{c}{p}_{00}=\frac{p_{10}}{2\gamma p_{12}+\gamma p_{11}+p_{10}}\\ p_{01}=\frac{\gamma p_{11}}{2\gamma p_{12}+\gamma p_{11}+p_{10}}\\ p_{02}=\frac{2\gamma p_{12}}{2\gamma p_{12}+\gamma p_{11}+p_{10}}\end{array}

Dominant

\begin{array}{l}{p}_{00}=\frac{p_{10}}{\gamma p_{12}+\gamma p_{11}+p_{10}}\hfill \\ p_{01}=\frac{\gamma p_{11}}{\gamma p_{12}+\gamma p_{11}+p_{10}}\hfill \\ p_{02}=\frac{\gamma p_{12}}{\gamma p_{12}+\gamma p_{11}+p_{10}}\hfill \end{array}

Recessive

\begin{array}{c}{p}_{00}=\frac{p_{10}}{\gamma p_{12}+p_{11}+p_{10}}\\ p_{01}=\frac{p_{11}}{\gamma p_{12}+p_{11}+p_{10}}\\ p_{02}=\frac{\gamma p_{12}}{\gamma p_{12}+p_{11}+p_{10}}\end{array}

After the power is computed for each SNP, I take the overall power to be the average power over all the SNPs. In taking the average power over all SNPs, I give less weight to the tag SNPs since they are over-represented in the set of SNPs being analyzed. Assume that of the S SNPs under consideration (for which we have linkage disequilibrium [LD] data from, for instance, the HapMap project), M are tags that will be genotyped on the chip and S-M are not tags. Further assume that there are T common SNPs in total in this population, which includes both those S SNPs for which we have LD data and SNPs for which we do not know their LD with surrounding SNPs. Let 1-β _i be the power for SNP i where i ranges from 1 to S and SNP i is a tag SNP when i ≤ M and a non-tag otherwise. Then, the overall power is given by:

Power={\displaystyle \sum _{i=1}^M\left[\frac{1}{T}(1-{\beta }_i)\right]}+{\displaystyle \sum _{i=M+1}^S\left[\frac{T-M}{T(S-M)}(1-{\beta }_i)\right]}

In this manner, the tag SNPs are only considered representative of themselves, while the non-tag SNPs for which we have LD data are considered representative of all common non-tag SNPs. For these calculations, I use T = 2 × 10⁷.

Implementation

A computer program to implement these calculations was written in C. The source code is available upon request from the author.