In this study, we examine statistical methods used to perform genome wide association studies (GWAS). GWAS usually apply univariate statistical tests to each gene marker or single nucleotide polymorphism (SNP) as an initial step. This SNP based test is statistically straightforward and the testing is done with standard methods (e.g. χ² tests, regression) that have been studied outside of the GWAS context. A recent paper by Kuo and Feingold (2010) described the most commonly used methods, and the authors note that a compound procedure which combined two or more statistical tests is used.

The literature contains a number of papers that make statistical power comparisons among subsets of these methods, including Sasieni (1997) and Freidlin et al. (2002), but the question of which method is best suited to univariate scanning in a GWAS remains an open issue. The choice of method depends on the match between the true genetic model underpinning the association and the type of model assumed by the method.

We used a multiple test procedure that combines the most promising of the methods identified in the literature and apply them to a set of synthetic marker data with known properties. Our goal is to identify marker properties that can be linked to optimal methods (with reference to statistical power) for predicting associations in GWAS. We know from prior studies that the statistical procedure a researcher chooses influences GWAS prediction accuracy, and that there are specific properties of the underlying markers that determine the optimization of the procedure choice, see Kuo and Feingold (2010). We included the important properties that influence the association prediction accuracy into our synthetic marker data via a Monte Carlo simulation process, and we link the properties to the influencing marker to study their individual and collective contributions to association prediction. A synthetic marker dataset allows us to assess the performance of different statistical methods in a GWAS context. We apply a number of statistical methods to the simulated data and use their statistical power profiles to evaluate the performance of the methods. We also quantify the relationships between locus traits and prediction accuracy.

Our study identifies a number of these properties and quantifies the loss in power if non-optimal methods are used. Similar results have been reported in earlier studies by Sasieni (1997) and Freidlin et al. (2002). These studies reinforce the view that the major influence on prediction accuracy is the gene model of the locus associated with the diagnosis.

We are particularly interested in assessing the consequences of applying a statistical method that assumes an inherent additive mode of inheritance (MOI) property to SNP data that is non-additive. Our motivation for this is twofold. First, the additive MOI model is commonly employed in GWAS, and second, the answer to the question: "what statistical methods should be used to conduct GWA type studies?" does not have a definitive answer. The best method typically depends on what MOI gene model has been associated with the associated diagnosis.

Our results show the major factors that influence association predictions. They also indicate that a strategy based on predicting associations using multiple statistical methods can be more accurate (much more accurate if the governing marker is recessive) than those that assume a single (additive) mode. The multiple test procedure proposed here uses a combination of the recessive, additive and dominant MOI-optimal statistical methods, all of which are derived from the well known Cochran-Armitage test. We also examined different procedures for combining the tests.

We examine the accuracy of association detection by generating synthetic data with properties that are known to influence statistical power. We used a Monte Carlo method to generate the data from a set of random variables described below. The main purpose of the synthetic data is to act as a "truth set" to assess the performance of commonly used statistical methods used in a GWAS context.

Generating the data

Our method for generating the synthetic marker data is derived from a study by Iles (2002) and is based on Mendelian concepts of inheritance. We include autosomal dominant and autosomal recessive patterns that are single gene inheritance patterns. We also incorporate additive and multiplicative inheritance patterns to represent the actions of multifactorial inheritance patterns.

A description of the data generating process begins from the notion of disease penetrance. Penetrance in genetics is the proportion of individuals carrying a particular variation of a gene (allele or genotype) that also express an associated trait. We designate A as the risk allele, and a as the allele without risk. By using the relationships between penetrance and relative risk as defined in Table 1, for different MOI categories, generating the synthetic gene dataset was straightforward. Specifically the steps were:

Table 1: Relative Risk assumptions by MOI, see Iles (2002).

1. Preload the details that define the factor combinations per MOI category (we used 864 combinations for this study). The factors are:

1.1 n_j = the target number of cases and controls in a given experiment (100, 250, 500,1000, 2000,4000, 8000, 9500),

1.2 dp_j = the disease penetrance, (.3, .4, .5),

1.3 ErrP_j = the misclassification error rate contained in the phenotype data, (0, 2, 5%),

1.4 ErrG_j = the misclassification error rate contained in the genotype data, (0, 2, 5%),

1.5 Φ_j = the relative risk, (1.0, 1.15, 1.3, 1.45),

2. Draw a genotype distribution (at random from the master set of genotype distributions obtained from real distribution data (i.e., the Schymick et al. (2007) data). At this stage, Chan et al. (2009) recommends that a minor allele frequency (MAF) threshold not be applied. They argue that filtering MAFs out of the process because of low frequencies or to maintain Hardy-Weinberg equilibrium (HWE) deviation has little effect on the overall false positive rate and in some cases, filtering MAF only serves to exclude SNPs. The effect of this step is to select a specific genotype distribution (at random) from the master distribution.

3. Use Table 1 to assign a case (1) or a control (0) based on the selected on genetic relative risk (GRR), penetrance (P) and MOI factors. This step converts the GRR ratio value into the probability that the case occurs for the MOI gene model of interest. This process can be represented by the following logic that was derived from Iles (2002):

Major Homozygote (aa): Assume that the aa genotype is selected. The probability of a case given this selection is equal to the disease penetrance P, or. ψ_aa = P.

Minor Homozygote (AA) - liability increasing allele: Assume the aa genotype is selected. GRR can be expressed as a ratio of two probabilities: the probability of a case for a minor homozygote divided by the probability of a case for a major homozygote, i.e.,

ψ_AA = Prob(case/AA) / Prob(case/aa) = x/P. (1)

From (1) the probability of a case given the minor genotype = x = ψ_AA * P (2)

where ψ_AA = one of the assigned risk factors and P is one of the assigned penetrance factors.

Heterozygote (aA): Assume the aA genotype is selected. By the same argument, the phenotype risk given a heterozygote is:

ψaA = Prob(case/aA) / Prob(case/aa) = y/P, (3)

By the same arguments, the risk of a case given the heterozygote = y = ψaA * P (4)

where ψaA = one of the assigned risk factors and P is one of the assigned penetrance factors.

4. Using the estimate of x and y, assign a case or control at random using the four different MOI models in conjunction with equations (2) and (4) and Table 1 below. We assigned cases in proportion to x (y) and controls in proportion to 1- x (1-y) for the minor homozygote (heterozygote) genotypes respectively. For the MOI models that assume an elevated risk from the minor and the hetero genotypes, we would expect a higher proportion of cases to be more easily identified via the statistical procedures. The specification of risk depends on specific and unknown disease mechanisms. A relative risk of 1.7 is considered strong and is associated with positive replication, see Sladek et al. (2007), and a risk of 1.3 is considered by Ziegler et al. (2008) to be a realistic assumption for complex diseases. We limited our focus to a relative risk range of 1.15 to 1.45 and were particularly interested in scenarios with low relative risk.

5. Individuals are either assigned as cases or controls according to the probabilities given in Table 1.

6. Continue with the above process until n1 cases and n2 controls are generated (note in this example n1 = n2, but that can be tailored to a specific set of n1 - n2 targets).

7. Apply a set of statistical methods to predict associations and record the results.

8. Generate 1,000 replicate experiments for each set of 3,456 factor combinations.

The simulated dataset that we generated has the following characteristics:

• The proportion of cases (controls) that are major homozygotes = 50.3 (63.0)%.

• The proportion of cases (controls) that are heterozygotes = 39.2 (31.3)%.

• The proportion of cases (controls) that are minor homozygotes = 10.5 (5.7)%.

• With MOI distribution:

ο Recessive = 25%,

ο Dominant = 25%,

ο Additive = 25%, and

ο Multiplicative = 25%.

We acknowledge that the distribution of MOI traits above is not representative of how inheritance traits are distributed in humans. The Online Mendelian Inheritance in Man (OMIM, http://www.ncbi. nlm.nih.gov/omim) provides the best source of information on the MOI distribution (Table 2). However, OMIM is disproportionally populated by genes linked to single Mendelian disorders. Therefore, genes associated with multifactorial disorders are under-represented in OMIM. Because polygene influences are assumed to be a major source of additive and multiplicative SNP behavior, the distribution in Table 2 is likely biased. Accordingly, we populated SNPs in our data with equal MOI representation and acknowledge that it is not representative of the true distribution.

Table 2: Distribution of Genes in OMIM by MOI.

The three "optimal" MOI specific methods are the three variations of the Cochran-Armitage (CA) trend test described in Zheng and Gastwirth (2006). We also included a fourth "individual" method, the 2df genotype test, which is a commonly used method. Using the notation in Table 3 to define the 2x3 table of case-control counts stratified by genotype , a test statistic (T²(x)) for the three variations of the CA trend methods is defined as:

T²(x) = n [Σ_0,1,2 {x_i (s r_i - r s_i)}]² / [r s (Σ_0,1,2 n {x_i x_i n_i } - {Σ_0,1,2 (x_i n_i)² })]. (5)

Table 3: Terms Defined in Equations (5) and (6).

The values represented in equations (5) and (6) are shown in Table 3 below and the value of x_i defines the specific test x0 = 0, x2 = 1 and x1 = {0 - recessive, .5 - additive, 1 - dominant).

Under the null hypothesis of no association, T²(x) has an asymptotic χ² distribution with 1 degree of freedom.

As an alternative to equation (5), it is also possible to use a normally distributed test statistic per Li et al. (2008):

N(x) = n^1/2 [Σ_0,1,2 {x_i (s r_i - r s_i)}] / [r s (n Σ_0,1,2 {x_i x_i n_i} - [Σ_0,1,2 {x_i n_i}²]^1/2.

Under the null hypothesis of no association, N(x) has an asymptotic normal distribution N(0,1), which suggests a one-sided test because the synthetic data assumes that the minor allele conveys the risk of phenotype.

We use eight sample size assumptions with equal numbers of cases and controls to perform our analysis, with N defined as the number of cases = 100, 250, 500, 1,000, 2,000, 4,000, 8,000, or 9,500. We estimate statistical power by statistical method and N using a significance threshold of α = 10^-7. In a GWA study, researchers usually perform a single marker analysis as a starting point to identify SNPs for additional and more comprehensive analysis. This initial pass creates a large number of statistical tests as well as a high potential for false-positive predictions, which has caused researchers to perceive the need for a very low threshold. Accordingly, recent studies have used type–I threshold levels on the order of 10^-7 as in Iles (2002), Ziegler et al. (2008) and Van Es et al. (2008).

The multi-test statistical methods we use in our comparisons are:

1. The Bonferroni (BON) method, shown in Holm (1979). A simple form of the Bonferroni correction results when using n methods to test for an association outcome. The correction involves dividing the alpha level by n. For example, if the association of a given SNP involves using three different statistical methods, the corrected alpha level (α) would be α/3. This would ensure that the overall chance of making a Type I error is still less than α.

2. A MAX method from Li et al. (2008) that departs from the Bonferroni method. Bonferroni assumes that the individual tests are mutually independent, while Li et a l. (2008) assumes that the individual tests are correlated and incorporates an approximation of the joint distributions into the method.

Results - statistical method assessment

Table 4 presents power estimates by statistical methods and sample size, and is based on a fixed alpha threshold (α = 10^-7). All tests are one-sided and the tests included in this table are:

Table 4: Power Results by Statistical Method and Number of Cases: Additive MOI Data.

• The additive χ² version of the CA (CA-A) test, which was the best method for both additive and multiplicative gene models but not particularly effective when applied to recessive MOI data.

• The Bonferroni test (BON), which uses (combines) the χ² version of the recessive, additive, and dominant MOI specific tests (CA-R, CA-A, CA-D) and improves on the test performance of the three individual tests when the MOI gene model is not known.

• The MAX test due to Li et al. (2008) which uses (combines) the Normal version of the CA-R, CA-A, CA-D tests to improve on the test performance of the three individual tests if the MOI gene model is not known.

• The dominant χ² version of the CA test (CA-D), which was the best method for the dominant gene models.

• The recessive χ² version of the CA test (CA-R), which was the best method for recessive gene models but the least effective when applied to non-recessive MOI data.

• The 2df genotype test (2df-g), which is never an optimal method for any of the scenarios we examined; in every scenario, a more powerful alternative can be identified (see Table 4).

Our results indicate that the best method in terms of statistical power is CA-A, but that little is lost if the MAX method is used instead. Similarly, the results in Tables 5, 6 and 7 indicate that the best method in terms of statistical power for identifying dominant MOI loci is CA-D and CA-R for recessive MOI loci. For multiplicative MOI loci, the best method is CA-A. In all four scenarios, little is lost if the MAX (or the BON) methods are used as replacements.

Table 5: Power Results by Statistical Method: Dominant MOI Gene Data.

Table 6: Power Results by Statistical Method: Recessive MOI Gene Data.

Table 7: Power results by Statistical Method: Multiplicative MOI Gene Data.

However, if we use the CA-A method, which is advocated by many as an initial GWAS pass and the locus in question is recessive or dominant, a power loss will occur. Table 8 indicates that while the CA-A method is the optimal choice (by as much as 2%), if the MOI of the locus is additive or multiplicative, there is a risk of as much as a 4% power loss if the locus MOI is dominant and as much as 15% loss if the MOI of the locus is recessive (see Table 8).

Table 8: Power results by CA-A and MAX Methods, for different MOI gene models.

If we knew the distribution of the MOI property, we could assess the overall risk of using an additive method such as CA-A for GWAS. However, without a reliable estimate, one should exercise caution and apply a procedure that limits the risk of making the wrong assessment of the MOI inherent to the locus inducing diagnosis.

In the literature, many statistical methods that have been used to perform GWAS assume a MOI specific hypothesis. Our results confirm the work of many others (see Iles, 2002): that is, in the context of a single marker scenario, the best method for predicting associations in recessive SNPs was the CA-R method; the best method for dominant MOI SNPs was the CA-D method; and the best method for additive and multiplicative SNPs was the CA-A method. We also show that the 2df genotype method used in many studies (for example, see Schymick et al. (2007)) is never an optimal method because there are always other methods that provide greater statistical power. This statement holds regardless of whether the MOI is known a priori or not. We also show that in the context of a general method to use in the initial GWAS pass, researchers may encounter adverse consequences if, for example, the MOI of the operating locus is not consistent with the assumption employed by the statistical method used. Thus, using 2df does not appear to be appropriate for GWA studies in any circumstances.

Consequently, we examined the possibility of employing an alternative procedure that incorporates the three core tests defined above into two multi test procedures: the Bonferroni procedure and the MAX test procedure due to Li et al. (2008). These procedures are opposites in many respects, in that they assume different underlying distributions of the test statistics. The MAX method assumes that the three separate tests have dependencies that can be accounted for, whereas the Bonferroni method assumes that they are mutually independent. We note that despite these differences, the two methods produce similar power profiles.

We generated our results using 1,000 replicates per parameter combination. Our standard error estimate of power varies from .262 to .315. Consequently, our 95 percent confidence interval around the mean will be approximately plus or minus .019. While we recognize that a larger number of replicates will improve power precision, we believe that our conclusion will remain as stated.

In summary, researchers should consider a multi-test procedure combining the results of individual MOI-based core tests as a possible statistical method for conducting GWAS rather than a 2df test. Combining individual methods and comparing the individual and combined results may help identify the MOI character of the gene. The actual process of combining the core tests into a single operational test can be done in a number of ways, all of which produce very similar statistical power profiles.

The authors wish to thank Craig R. Hollingsworth of RTI for technical writing and editing assistance.