 Methodology article
 Open Access
 Published:
Calculation of exact pvalues when SNPs are tested using multiple genetic models
BMC Geneticsvolume 15, Article number: 75 (2014)
Abstract
Background
Several methods have been proposed to account for multiple comparisons in genetic association studies. However, investigators typically test each of the SNPs using multiple genetic models. Association testing using the CochranArmitage test for trend assuming an additive, dominant, or recessive genetic model, is commonly performed. Thus, each SNP is tested three times. Some investigators report the smallest pvalue obtained from the three tests corresponding to the three genetic models, but such an approach inherently leads to inflated type 1 errors. Because of the small number of tests (three) and high correlation (functional dependence) among these tests, the procedures available for accounting for multiple tests are either too conservative or fail to meet the underlying assumptions (e.g., asymptotic multivariate normality or independence among the tests).
Results
We propose a method to calculate the exact pvalue for each SNP using different genetic models. We performed simulations, which demonstrated the control of type 1 error and power gains using the proposed approach. We applied the proposed method to compute pvalue for a polymorphism eNOS 786T>C which was shown to be associated with breast cancer risk.
Conclusions
Our findings indicate that the proposed method should be used to maximize power and control type 1 errors when analyzing genetic data using additive, dominant, and recessive models.
Background
Genomewide association studies (GWAS) and candidate gene association studies are commonly performed to test the association of genetic variants with a particular phenotype. Typically, hundreds of thousands of singlenucleotide polymorphisms (SNPs) are tested for association in these studies. Associations between the SNPs and the phenotypes are determined on the basis of differences in allele frequencies between cases and controls [1]. Several statistical methods have been proposed to control the familywise error rate (FWER) for multiple comparison testing.
A simple approximation can be used to obtain a FWER of α by utilizing the Bonferroni adjustment [2] of ${\mathit{\alpha}}^{*}=\frac{\mathit{\alpha}}{\mathit{n}}$ and using α* as the threshold for significance for each test. Bonferroni adjustment tends to be conservative when the tests are correlated. In genetic association studies, the SNPs being tested are typically in linkage disequilibrium (LD), which leads to correlation among the tests. An alternative approximation to the Bonferroni adjustment is Sidak’s correction [3, 4], ${\mathit{\alpha}}^{*}=1{\left(1\mathit{\alpha}\right)}^{\frac{1}{\mathit{n}}}$ which assumes independence among tests. Conneely and Boehnke [5] proposed a correction that does not assume independence among tests but assumes joint multivariate normality of all test statistics. Other methods to control the FWER include using the false discovery rate (FDR) [6, 7].
In genetic association studies, three genetic modelsadditive, dominant, and recessiveare generally used to test each SNP using the CochranArmitage (CA) trend test [8–12]. In association studies the true underlying genetic model is unknown. Some investigators report the smallest pvalue obtained from the three tests corresponding to the three genetic models. However, such a procedure inherently leads to an inflated type 1 error rate. Also, FDRbased methods to control FWER are not applicable in this situation because the hypotheses are highly correlated, as the same SNP is tested using different genetic models.
Thus, there is a need to correct for multiple comparisons corresponding to the three genetic tests performed for testing the association of a single SNP. These three tests are not only correlated but also functionally dependent. The standard methods for correcting for multiple testing referred to above are either too conservative or fail to meet the assumptions underlying these methods (e.g., asymptotic multivariate normality, independence among tests). Several approaches have been proposed to account specifically for the multiple comparisons of these three genetic models [13–15]. However, these approaches assume asymptotic trivariate normality for the additive, dominant and recessive test statistics. While this is a reasonable approximation to correct for multiple comparisons, preliminary investigations regarding the joint distribution of the three test statistics revealed the following insights: 1) the joint distribution of the test statistics is discrete and the grids at which the probability mass function is positive is few and far between; 2) The distribution is highly multimodal in most of the situations, particularly, when the number of cases and controls are different and unimodal only in a handful of situations (e.g. when the number of cases and controls are equal). Therefore, we propose a method to compute the exact joint distribution of the three CA trend tests corresponding to the additive, dominant, and recessive genetic models. We used this joint distribution to compute the exact pvalue for testing each SNP using the different genetic models. We performed simulations to demonstrate control of type 1 errors and power gains using the proposed approach. Finally, we applied the proposed approach to assess the significance of the association between a promoter polymorphism, eNOS786T>C and breast cancer risk.
Methods
Consider a diallelic SNP locus. The minor (deleterious) allele is labeled as a, and the major (normal) allele is labeled as A. The deleterious allele a is assumed to affect a phenotype Z, which takes the values of 0 or 1: Z = 1 indicates cases (affected) and Z = 0 indicates controls (unaffected). The observed genotype data for the SNP is one of three genotypes (A, A), (A, a), or (a, a). Let R_{ X } denote the number of cases and R_{ Y } denote the number of controls, with R_{ X } + R_{ Y } = N. Let X_{1}, X_{2}, X_{3} and Y_{1}, Y_{2}, Y_{3} be the number of individuals with genotypes AA, Aa, and aa in cases and controls, respectively. The data can be formulated in a 2 × 3 contingency table, as shown in Table 1. Let p_{1}, p_{2}, p_{3} be the frequencies of genotypes, AA, Aa and aa in cases and q_{1}, q_{2}, q_{3} be the frequencies of these three genotypes in controls. The values of p_{ i }, q_{ i }, i=1,2,3 can be estimated from the data as $\phantom{\rule{0.25em}{0ex}}{\mathit{p}}_{\mathit{i}}=\frac{{\mathit{X}}_{\mathit{i}}}{{\mathit{R}}_{\mathit{X}}}$ and $\phantom{\rule{0.25em}{0ex}}{\mathit{q}}_{\mathit{i}}=\frac{{\mathit{Y}}_{\mathit{i}}}{{\mathit{R}}_{\mathit{Y}}}\phantom{\rule{0.25em}{0ex}}$.
There have been many approaches in the literature for testing the association between a SNP and disease status. The CA test for trend [8] is generally the most popular and is available in most genetic analysis software packages, such as PLINK [16]. The test statistic for the CA test is as follows:
where the weight, t_{ i }, is chosen on the basis of the genetic model considered: additive, dominant, or recessive. The additive model assumes the deleterious effect is linearly related to the number of deleterious alleles. The dominant model assumes the deleterious effect is related to the presence of the deleterious allele. And the recessive model assumes the deleterious effect is related to the presence of both the deleterious alleles. The weights t = [t_{1}, t_{2}, t_{3}] corresponding to each of the models are as follows: additive model: t = [0, 1, 2], dominant model: t = [0, 1, 1], and recessive model: t = [0, 0, 1] for genotypes AA, Aa, and aa, respectively. Let the three test statistics corresponding to the additive, dominant, and recessive models be T_{1}, T_{2}, and T_{3}, respectively.
The joint distribution
Each test statistic, T_{1}, T_{2} and T_{3}, has an asymptotically normal univariate distribution. Therefore, the pvalues for each of these tests can be obtained from their asymptotic distributions. However, reporting the smallest pvalue obtained from testing T_{1}, T_{2} and T_{3,} individually leads to an inflated type 1 error rate. If the exact joint distribution of the three tests is known, one can compute the exact pvalue for the SNP that will account for the multiple correlated tests. We proceed to derive the joint distribution of the three test statistics, T_{1} = (R_{ Y }X_{2}−R_{ X }Y_{2}) + 2(R_{ Y }X_{3} − R_{ X }Y_{3}), and T_{2} = (R_{ Y }X_{2} − R_{ X }Y_{2}) + (R_{ Y }X_{3} − R_{ X }Y_{3}), and T_{3} = (R_{ Y }X_{3} − R_{ X }Y_{3}). As T_{3} = T_{1} − T_{2}, we only need to derive the joint distribution of T_{1} and T_{2}. It is reasonable to assume that the three genotype counts in cases (X_{1}, X_{2}, X_{3}) and the three genotype counts in controls (Y_{1}, Y_{2}, Y_{3}) follow a multinomial distribution, with probabilities (p_{1}, p_{2}, p_{3}) and (q_{1}, q_{2}, q_{3}) respectively. Let $\mathit{T}=\left(\begin{array}{c}\hfill {\mathit{T}}_{1}\hfill \\ \hfill {\mathit{T}}_{2}\hfill \end{array}\right)$, $\mathit{X}=\left(\begin{array}{c}\hfill {\mathit{X}}_{2}\hfill \\ \hfill {\mathit{X}}_{3}\hfill \end{array}\right)$ and $\mathit{Y}=\left(\begin{array}{c}\hfill {\mathit{Y}}_{2}\hfill \\ \hfill {\mathit{Y}}_{3}\hfill \end{array}\right)$. The test statistics can be written as T = AX + BY, where $\mathit{A}=\left[\begin{array}{cc}\hfill {\mathit{R}}_{\mathit{Y}}\hfill & \hfill 2{\mathit{R}}_{\mathit{Y}}\hfill \\ \hfill {\mathit{R}}_{\mathit{Y}}\hfill & \hfill {\mathit{R}}_{\mathit{Y}}\hfill \end{array}\right]$ and $\mathit{B}=\left[\begin{array}{cc}\hfill {\mathit{R}}_{\mathit{X}}\hfill & \hfill 2{\mathit{R}}_{\mathit{X}}\hfill \\ \hfill {\mathit{R}}_{\mathit{X}}\hfill & \hfill {\mathit{R}}_{\mathit{X}}\hfill \end{array}\right]$. Then the joint probability mass function (pmf) of T_{1}, T_{2} is given by
where f_{ x }, f_{ y } are trinomial probability mass functions and h(X, T) = B^{−1}T − B^{−1}AX. The derivation of the joint pmf of T_{1}, T_{2} is detailed in the Appendix. The pvalue corresponding to the test statistic (t_{1}, t_{2}) can be computed by summing up the probabilities of the test statistics that are equally or less probable than the observed test statistic, which can be written as
The computation of the pvalue using the above formula is nontrivial; however, there are a variety of computational optimizations and parallels to Fisher’s exact test that can be used to drastically reduce the computational complexity (see details in the Appendix). Briefly, the CA trend test statistics form a system of constrained linear Diophantine equations. The computational optimizations presented in the Appendix are based on exploiting the properties of the linear Diophantine equations with trinomial constraints. The solution space of these equations corresponds to the discrete space of nonzero probabilities for the joint pmf. This discrete space has a pattern of overlapping triangles that can be enumerated based on R_{ X } and R_{ Y } counts (See Figures 1, 2, 3 and 4). To reduce the number of computations in the discrete space we first transformed the test statistics to be symmetric. The pattern of overlapping triangles depends on three different scenarios based on the greatest common divisor (GCD) of R_{ X } and R_{ Y }: 1. GCD(R_{ X }, R_{ Y }) = 1, 2. GCD(R_{ X }, R_{ Y }) = R_{ X } = R_{ Y } and 3. 1 < GCD(R_{ X }, R_{ Y }) < min(R_{ X }, R_{ Y }). In scenario 1 the triangles do not overlap, therefore the pvalue can be evaluated most efficiently (Figures 1 and 2). In scenario 2 most of the triangles overlap and the discrete space of nonzero probabilities is sparse (Figure 4). In this scenario, we proposed an algorithm to exploit this aspect to calculate the exact pvalue more efficiently. Scenario 3 is the most general case which uses the general optimizations of symmetricity and the triangle pattern (Figure 3). The algorithms to compute the exact pvalues for each of the scenarios are detailed in the Appendix.
Simulations
We performed simulations to evaluate the performance of the proposed method and compared our approach with standard approaches used in the literature. All the simulation results were based on 1000 replicate data sets. Each replicate dataset comprised 1000 cases and 1000 controls. The disease status for each data set was obtained using the logistic regression model logit(P(Z = 1)) = β_{0} + β_{1}X, where X is the indicator for genotype, Z is the disease status, β_{0} is the intercept, and β_{1} is the log odds ratio for the SNP. The genotype data for a SNP were simulated using a minor allele frequency (MAF) of 40% for the null hypothesis and two MAFs of 40% and 20% for the power comparisons. For the type 1 error comparisons, we simulated 1000 replicate datasets from the null hypothesis (i.e., the SNP was not associated with disease status), with β_{0} = − 2.5 and β_{1} = log (1). For the power comparisons, we simulated 1000 replicate datasets for 40% and 20% MAFs from the alternate hypothesis (i.e., the SNP was associated with disease status) for each of the three scenarios: (1) additive model with odds ratio of 1.2, (2) dominant model with odds ratio of 1.3, and (3) recessive model with odds ratio of 1.3. The methods we compared were as follows: performing only additive analyses (additiveonly), performing only dominant analyses (dominantonly), performing only recessive analyses (recessiveonly), using the pvalue based on reporting the smallest pvalue of the three genetic models (minp), using the Bonferroni correction approach, and using the proposed exact pvalue method.
Results
The type 1 errors based on 1000 replicates from the null hypothesis are shown in Table 2. Analyses based on additiveonly, dominantonly, and recessiveonly models gave empirical type 1 errors of 0.044, 0.045, and 0.056, respectively, at the 0.05 level of significance. As expected, these models provided good control of type 1 errors because only one genetic model was tested in these analyses. The Bonferroni approach also had a wellcontrolled, but conservative, type 1 error (0.030 at the 0.05 level of significance). The minp had a type 1 error of 0.105 at the 0.05 level of significance, which was very liberal and confirmed that the minimum pvalue of the three genetic models is not a valid test. Finally, our proposed approach provided good control of the type 1 error (0.047 at the 0.05 level of significance).
The power comparisons based on 1000 replicates for the SNP data simulated using 40% and 20% MAFs for the three scenarios when the data were simulated using the additive, dominant, and recessive models, respectively, are shown in Table 3. The top and bottom panels of Table 3 depict the results for 40% and 20% MAFs, respectively. The minp model was excluded from the comparison because of its inflated type 1 error. When the data were simulated using the additive genetic model (column 3, Table 3), and were analyzed using only the additive model, it had the highest powers (0.816 and 0.656 for 40% and 20% MAFs, respectively). However, when the data were analyzed using only the dominant model, the powers were 0.676 and 0.603 for 40% and 20% MAFs, respectively. Also, when the data were analyzed using only the recessive model the powers were 0.588 and 0.306 for 40% and 20% MAFs, respectively. The powers for the additive only analysis were the highest as expected because the true simulation model in this scenario was additive. However, the true model of disease inheritance is generally unknown and one performs analyses using all three genetic models. In this scenario, the proposed exact pvalue method had powers of 0.743 and 0.584 for 40% and 20% MAFs, respectively, at the 0.05 level of significance, which were higher than the Bonferroni method which had powers of 0.721 and 0.556 for 40% and 20% MAFs, respectively. Overall, powers of the proposed method were lower than additive model (true simulation model) but higher than those of the dominantonly, recessiveonly, and Bonferroni correction approach.
When the data were simulated using the dominant model (column 4, Table 3), the additiveonly, dominantonly and recessiveonly analyses had powers of 0.660, 0.803, and 0.158, respectively, for 40% MAF and 0.774, 0.823, and 0.102, respectively for 20% MAF, at the 0.05 level of significance. Once again, as expected, the powers of the dominantonly analysis were the highest because the data were generated using the dominant model. The proposed exact pvalue method had powers of 0.726 and 0.782 for the 40% and 20% MAFs, respectively, which were higher than the Bonferroni method which had powers of 0.671 and 0.715 for the 40% and 20% MAFs, respectively. When the data were simulated using the recessive model (column 5, Table 3), the additiveonly, dominantonly and recessiveonly analyses had powers of 0.410, 0.116, and 0.589, respectively, for 40% MAF and 0.116, 0.061, and 0.249, respectively, for 20% MAF. The proposed exact pvalue method had powers of 0.517 and 0.197 for the 40% and 20% MAFs, respectively, which were higher than the Bonferroni method (0.452 and 0.168 for 40% and 20% MAFs, respectively).
We applied the proposed approach to assess the significance of the association between the promoter polymorphism eNOS 786T>C and sporadic breast cancer risk in nonHispanic white women younger than 55 years from a breast cancer study performed by [17]. The study discovered that eNOS 786T>C was statistically significant for breast cancer (p=0.017) and included 421 breast cancer cases and 423 cancer free controls. The first panel in Table 4 depicts the genotype counts for TT, CT and CC genotypes in cases and controls for the eNOS 786T>C. The second panel in Table 4 reports the pvalues for the eNOS 786T>C computed using the 5 different approaches: additiveonly, dominantonly, recessiveonly, Bonferroni and the proposed exact pvalue method. The additiveonly, dominantonly and recessiveonly approaches had pvalues of 0.0045, 0.0148 and 0.0313, respectively, and the Bonferroni adjusted pvalue was 0.0135. For this SNP, the pvalue computed using the proposed exact pvalue method was 0.0021, which was more significant than the smallest of the three pvalues obtained using the additive, dominant, and recessiveonly analyses (Table 4).
Discussion
In this paper, we proposed a method to calculate the exact pvalue for testing a single SNP using multiple genetic models. We recommend using the proposed method to maximize power and control type 1 errors when analyzing genetic data using additive, dominant, and recessive models. The proposed method is robust to model misspecifications and different SNP minor allele frequencies. Furthermore, similar to the computation of Fisher’s exact pvalue, the proposed approach does not depend on asymptotic distributions.
In our simulation study, where replicate datasets were simulated using the null hypothesis, we found that the proposed method had wellcontrolled type 1 error probabilities. In contrast, the method of reporting the smallest pvalue of the three genetic models tested had the highest falsepositive rate and was found to be invalid. And, as expected, the type 1 error of the Bonferroni correction approach was well controlled but conservative, which typically led to a loss in power for identifying genetic variants.
We also simulated replicate datasets under an alternative hypothesis using the different genetic models: additive, dominant, and recessive. In these simulations, we observed that no single method: additiveonly, dominantonly, or recessiveonly, had higher power in all three scenarios. Each of these methods had higher power only when the model used to analyze the data was the same as the true model used to generate the data. However, because the true mode of disease inheritance is usually unknown, analyses using all three genetic models are necessary. In general, the Bonferroni correction approach led to higher power than using a model that did not correspond to the true model. The proposed exact pvalue method was an improvement over the Bonferroni method. The conservativeness of the Bonferroni method may be due to its inability to account for the functional dependence between the three test statistics. In contrast, our proposed approach accounts for this functional dependence by computing pvalues from the joint probability mass function. Finally, we analyzed breast cancer study data in which the polymorphism eNOS 786T>C, was found to be significant [17].
The computation time needed to obtain the exact pvalue is substantial. The problem is very closely related to Fisher’s exact test, and there are many patterns inherent in the structure of the problem that could be exploited to calculate the pvalues more efficiently. In the Appendix, we present several novel optimization techniques to efficiently compute the test statistics in a reasonable time (e.g., approximately 15 min for a 1000 cases and 1000 controls dataset). The software to compute exact pvalues is available at http://odin.mdacc.tmc.edu/~rtalluri/index.html.
Conclusions
In genetic association studies, three genetic modelsadditive, dominant, and recessiveare generally used to test each SNP using the CochranArmitage trend test. Reporting the minimum pvalue of the three genetic models leads to inflated type 1 errors. We proposed an approach to compute the exact pvalue when genomic data is analyzed using the three genetic models. The proposed approach leads to higher power while controlling the type 1 error.
Appendix
Optimization techniques for computing the exact pvalue
Recall that X_{1}, X_{2}, X_{3} and Y_{1}, Y_{2}, Y_{3} are the number of individuals with genotypes AA, Aa, and aa in cases and controls, respectively, with X_{1} + X_{2} + X_{3} = R_{ X } and Y_{1} + Y_{2} + Y_{3} = R_{ Y }. The three genotype counts in cases (X_{1}, X_{2}, X_{3}) and the three genotype counts in controls (Y_{1}, Y_{2}, Y_{3}) follow a multinomial distribution with probabilities (p_{1}, p_{2}, p_{3}) and (q_{1}, q_{2}, q_{3}), respectively. The probability mass function (pmf) of (X_{ 1 }, X_{ 2 }, X_{ 3 }) is ${\mathit{f}}_{\mathit{X}}\left(\mathit{X}\right)=\frac{{\mathit{R}}_{\mathit{X}}!}{{\mathit{X}}_{2}!{\mathit{X}}_{3}!\left({\mathit{R}}_{\mathit{X}}{\mathit{X}}_{2}{\mathit{X}}_{3}\right)!}{\mathit{p}}_{1}^{{\mathit{R}}_{\mathit{X}}{\mathit{X}}_{2}{\mathit{X}}_{3}}{\mathit{p}}_{2}^{{\mathit{X}}_{2}}{\mathit{p}}_{3}^{{\mathit{X}}_{3}}$ and the pmf of (Y_{1}, Y_{2}, Y_{3}) is ${\mathit{f}}_{\mathit{Y}}\left(\mathit{Y}\right)=\frac{{\mathit{R}}_{\mathit{Y}}!}{{\mathit{Y}}_{2}!{\mathit{Y}}_{3}!\left({\mathit{R}}_{\mathit{Y}}{\mathit{Y}}_{2}{\mathit{Y}}_{3}\right)!}{\mathit{q}}_{1}^{{\mathit{R}}_{\mathit{Y}}{\mathit{Y}}_{2}{\mathit{Y}}_{3}}{\mathit{q}}_{2}^{{\mathit{Y}}_{2}}{\mathit{q}}_{3}^{{\mathit{Y}}_{3}}$. The three test statistics corresponding to the additive, dominant, and recessive models are, T_{1} = (R_{ Y }X_{2} − R_{ X }Y_{2}) + 2(R_{ Y }X_{3} − R_{ X }Y_{3}) , T_{2} = (R_{ Y }X_{2} − R_{ X }Y_{2}) + (R_{ Y }X_{3} − R_{ X }Y_{3}), and T_{3} = (R_{ Y }X_{3} − R_{ X }Y_{3}) respectively. As T_{3} = T_{1} − T_{2}, we only need to derive the joint distribution of T_{1} and T_{2}. Let $\mathit{T}=\left(\begin{array}{l}{\mathit{T}}_{1}\\ {\mathit{T}}_{2}\end{array}\right)$, $\mathit{X}=\left(\begin{array}{l}{\mathit{X}}_{2}\\ {\mathit{X}}_{3}\end{array}\right)$, and $\mathit{Y}=\left(\begin{array}{l}{\mathit{Y}}_{2}\\ {\mathit{Y}}_{3}\end{array}\right)$. The test statistics can be written as T = AX + BY, where $\mathit{A}=\left[\begin{array}{l}{\mathit{R}}_{\mathit{Y}}\phantom{\rule{0.25em}{0ex}}2{\mathit{R}}_{\mathit{Y}}\\ {\mathit{R}}_{\mathit{Y}}\phantom{\rule{0.25em}{0ex}}{\mathit{R}}_{\mathit{Y}}\end{array}\right]$ and $\mathit{B}=\left[\begin{array}{l}{\mathit{R}}_{\mathit{X}}\phantom{\rule{0.25em}{0ex}}2{\mathit{R}}_{\mathit{X}}\\ {\mathit{R}}_{\mathit{X}}\phantom{\rule{0.25em}{0ex}}{\mathit{R}}_{\mathit{X}}\end{array}\right]$ . We proceed to derive the joint probability mass function of $\mathit{T}=\left(\begin{array}{l}{\mathit{T}}_{1}\\ {\mathit{T}}_{2}\end{array}\right)$.
Consider an ndimensional discrete random vector G with pmf f_{ G }(). Suppose we have a transformation from G → H. The pmf f_{ H }() of the transformed variables H can be expressed as follows: [18]
This can be extended to the case where the dimensions of G and H are different, i.e., the transformation from (X, Y) → T is a linear transformation of the form T = AX + BY. The pmf of T is given by
This can be simplified as:
Computing this pmf on all the possible values of (T_{1}, T_{2}) is prohibitively time consuming. Computational optimizations can be used to speed up the computations of the probability mass function. We list several optimization techniques below. The first optimization is to transform the pmf to be symmetric in (T_{1}, T_{2}), which reduces the computational burden by half. The original test statistics T_{1} and T_{2} are T_{1} = (R_{ Y }X_{2} − R_{ X }Y_{2}) + 2(R_{ Y }X_{3} − R_{ X }Y_{3}) and T_{2} = (R_{ Y }X_{2} − R_{ X }Y_{2}) + (R_{ Y }X_{3} − R_{ X }Y_{3}), respectively. The joint pmf of (T_{1}, T_{2}) is a onetoone function of the joint distribution of any two orthogonal linear combinations of T_{1} and T_{2}. So if we transform the test statistics T_{1} and T_{2} into
the resulting pmf of (Z_{1}, Z_{2}) is a onetoone function of the pmf of (T_{1}, T_{2}). Hence, the pvalue obtained will be the same when using (Z_{1}, Z_{2}) instead of (T_{1}, T_{2}). The resulting pmf of (Z_{1}, Z_{2}) can be derived using the same method as with (T_{1}, T_{2}).
The next computational optimization is to identify the values that can be taken by (Z_{1}, Z_{2}). The number of values (Z_{1}, Z_{2}) can take are finite and represented by the solution space of the equations
which depends on the values of R_{ X } and R_{ Y }. These equations are called linear Diophantine equations and have an infinite number of solutions [19]. But in our case we have multiple constraints on the equations, which reduce the solution space to a finite number of solutions. The constraints are

1.
X _{3}, Y _{3}, X _{2} and Y _{2} are integers

2.
X _{3}, Y _{3}, X _{2} and Y _{2} ≥ 0

3.
X _{3} + X _{2} ≤ R _{ X }

4.
Y _{3} + Y _{2} ≤ R _{ Y }
On the basis of these four constraints the solution space can be calculated. While the exact solution space could not be found, it follows a pattern that can be enumerated.
Figure 1 depicts the pmf of the scenario with R_{ X } = 19 and R_{ Y } = 2 where a pattern of six triangles can be visualized from the figure. Similarly, Figure 2 depicts the pmf of the scenario with R_{ X } = 20 and R_{ Y } = 3, where a pattern of ten triangles can be visualized from the picture. This trend can be generalized for all values of R_{ X } and R_{ Y }.
Generalizing the above scenario, there are $\left[1+2+\cdot \cdot \cdot +\left({\mathit{R}}_{\mathit{Y}}+1\right)=\frac{\left({\mathit{R}}_{\mathit{Y}}+1\right)\left({\mathit{R}}_{\mathit{Y}}+2\right)}{2}\right]$ triangles for the solution space. In each triangle, there are $\left[1+2+\cdot \cdot \cdot +\left({\mathit{R}}_{\mathit{X}}+1\right)=\frac{\left({\mathit{R}}_{\mathit{X}}+1\right)\left({\mathit{R}}_{\mathit{X}}+2\right)}{2}\right]$ elements that correspond to all possible combinations of X_{3} + X_{2} ≤ R_{ X }. In each triangle, the values of Y_{3} and Y_{2} are constant and the $\frac{\left({\mathit{R}}_{\mathit{Y}}+1\right)\left({\mathit{R}}_{\mathit{Y}}+2\right)}{2}$ triangles correspond to all possible combinations of Y_{3} + Y_{2} ≤ R_{ Y }, which make up the whole solution space.
Another important fact is that these triangles may overlap, reducing the solution space, which is depicted in Figures 3 and 4. Figure 3 depicts the pmf of the scenario with R_{ X } = 10 and R_{ Y } = 2 where a pattern of six triangles can be visualized from the figure. The overlap of the triangles can be observed when compared to Figure 1. Figure 4 depicts the pmf of the scenario with R_{ X } = 5 and R_{ Y } = 5 where a pattern of 21 triangles can be visualized from the figure, where most of the triangles are overlapping one another. The additional computational burden is to determine where the solution space triangles overlap and how many triangles are overlapping at a particular location. This is a function of the greatest common divisor (GCD) of R_{ X } and R_{ Y }. If R_{ X } and R_{ Y } are coprime (GCD=1), only three triangles overlap at a single point (Z_{1} = 0, Z_{2} = 0) which requires no additional computation. When R_{ X } and R_{ Y } are not coprime, the triangles overlap at multiples of the GCD of R_{ X } and R_{ Y }. In this scenario, multiple values of X_{3}, Y_{3}, X_{2,} and Y_{2} contribute to the same (Z_{1}, Z_{2}).
In an ideal scenario, the total number of computations required to compute the pmf of (Z_{1}, Z_{2}) is $\frac{\left({\mathit{R}}_{\mathit{Y}}+1\right)\left({\mathit{R}}_{\mathit{Y}}+2\right)}{2}\frac{\left({\mathit{R}}_{\mathit{X}}+1\right)\left({\mathit{R}}_{\mathit{X}}+2\right)}{2}\approx \frac{{\mathit{R}}_{\mathit{X}}^{2}{\mathit{R}}_{\mathit{Y}}^{2}}{4}$, which can be computed in approximately 15 minutes for R_{ X } = 1000 and R_{ Y } = 1000 using a computer with a 3.4GHz processor and 8 GB of RAM. However, the amount of storage required for the solution space far exceeds the hardware capabilities available. In light of this limitation, computational optimizations should be employed to avoid storing the whole solution space. This limitation leads to three possible scenarios:

1.
GCD(R _{ X }, R _{ Y }) = 1

2.
GCD(R _{ X }, R _{ Y }) = R _{ X } = R _{ Y }

3.
GCD(R _{ X }, R _{ Y }) < min(R _{ X }, R _{ Y })
Scenario 1
When R_{ X } and R_{ Y } are coprime, the triangles only overlap at a single point (Z_{1} = 0, Z_{2} = 0); therefore, we can independently evaluate each of the possible values of the solution space. The pvalue is the probability of obtaining a test statistic at least as extreme as the one observed, so we evaluate the probabilities of each of the possible values of the test statistics one at a time. Hence, the pvalue is the sum of all the probabilities of test statistics that are lower than the probability of the observed test statistic. Using this procedure there is no need to store any data, which leads to faster computation of the pvalue from the joint distribution.
Scenario 2
When R_{ X } and R_{ Y } are equal, most of the triangles overlap with each other. But a pattern has been observed in this scenario, which is shown in Figure 5, where R_{ X } = 10 and R_{ Y } = 10. As seen in Figure 4, the solution space is very sparse and only requires computation of the colored cells. The possible solution space is spaced R_{ X } apart. So if we condense the possible solution space, the solution space is as shown in Figure 5. Figure 5 shows the number of triangles overlapping at each point in the solution space. Only half of the matrix needs to be computed, as the other half is symmetric. The algorithm to compute the pvalue is as follows.
Algorithm:

1.
Let R _{ X } = R _{ Y } = R. The solution space can then be constrained to a matrix with 2R + 1 rows and 2R + 1 columns. Let the center of the matrix correspond to the test statistic (Z _{1} = 0, Z _{2} = 0).2. Now, as we can see from Figure 5, we need to compute the colored cells in quadrants 3 and 4. In quadrant 3, the cells with the same number of overlapping triangles are placed diagonally, and in quadrant 4, they are placed horizontally and then vertically. We exploit the pattern that follows from the same number of triangles overlapping at a particular cell.

3.
For i = 1: R start at (Z _{1} = − (R − i), Z _{2} = − 1). Find the possible combinations of X _{3}, Y _{3}, X _{2} and Y _{2} that contribute to the cell corresponding to (Z _{1} = − (R − i), Z _{2} = − 1). Compute the probabilities for the cells along the diagonal path in quadrant 3, until Z _{1} = 0. Here X _{3} and X _{2} remain the same; hence, it is trivial to compute the probabilities for each cell.

4.
Then in quadrant 4, compute the probabilities for the cells along the horizontal path until Z _{1} = R − (i − 1); here X _{3} remains the same and X _{2new} = X _{2} + Z _{2}.

5.
Then continue vertically until Z _{2} = 0; here X _{3} and X _{2} remain the same.
This algorithm reduces the computational burden by computing the possible combinations of X_{3}, Y_{3}, X_{2} and Y_{2} that contribute to all the cells only R times, as opposed to computing once for each cell (approximately 4R^{2} times).
Scenario 3
This is the general scenario where GCD(R_{ X }, R_{ Y }) < min(R_{ X }, R_{ Y }). Several patterns that can be used to reduce the computational burden that could be applied for a particular GCD were found, but these could not be generalized to all the possible situations. We instead use a straightforward approach to determine the pvalue for each of the possible solutions for (Z_{1}, Z_{2}). The algorithm is as follows:

1.
For each possible (Z _{1}, Z _{2}) compute the triangles that contribute to this particular point.

2.
Add up the probabilities of each of the elements of these triangles to compute the pvalue of that particular (Z _{1}, Z _{2}).
References
 1.
Clarke GM, Anderson CA, Pettersson FH, Cardon LR, Morris AP, Zondervan KT: Basic statistical analysis in genetic case–control studies. Nat Protoc. 2011, 6 (2): 121133. 10.1038/nprot.2010.182.
 2.
Dunn OJ: Multiple comparisons among means. J Am Stat Assoc. 1961, 56 (293): 5264. 10.1080/01621459.1961.10482090.
 3.
Sidak Z: On multivariate normal probabilities of rectangles  their dependence on correlations. Ann Math Stat. 1968, 39 (5): 14251434.
 4.
Sidak Z: Probabilities of rectangles in multivariate student distributions  their dependence on correlations. Ann Math Stat. 1971, 42 (1): 169175. 10.1214/aoms/1177693504.
 5.
Conneely KN, Boehnke M: So many correlated tests, so little time! Rapid adjustment of P values for multiple correlated tests. Am J Hum Genet. 2007, 81 (6): 11581168. 10.1086/522036.
 6.
Benjamini Y, Hochberg Y: Controlling the false discovery rate  a practical and powerful approach to multiple testing. J Roy Stat Soc B Methods. 1995, 57 (1): 289300.
 7.
Benjamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. Ann Math Stat. 2001, 29 (4): 11651188. 10.1214/aos/1013699998.
 8.
Agresti A: Categorical Data Analysis. 2002, New York: John Wiley & Sons
 9.
Armitage P: Tests for linear trends in proportions and frequencies. Biometrics. 1955, 11 (3): 375386. 10.2307/3001775.
 10.
Barrett JH, Iles MM, Harland M, Taylor JC, Aitken JF, Andresen PA, Akslen LA, Armstrong BK, Avril MF, Azizi E, Bakker B, Bergman W, BianchiScarra G, Bressacde Paillerets B, Calista D, CannonAlbright LA, Corda E, Cust AE, Debniak T, Duffy D, Dunning AM, Easton DF, Friedman E, Galan P, Ghiorzo P, Giles GG, Hansson J, Hocevar M, Hoiom V, Hopper JL, et al: Genomewide association study identifies three new melanoma susceptibility loci. Nat Genet. 2011, 43 (11): 11081113. 10.1038/ng.959.
 11.
Cochran WG: Some methods for strengthening the common X2 tests. Biometrics. 1954, 10 (4): 417451. 10.2307/3001616.
 12.
Lewis CM, Knight J: Introduction to genetic association studies. Cold Spring Harb Protoc. 2012, 2012 (3): 297306.
 13.
Freidlin B, Zheng G, Li ZH, Gastwirth JL: Trend tests for case–control studies of genetic markers: power, sample size and robustness. Hum Hered. 2002, 53 (3): 146152. 10.1159/000064976.
 14.
Gonzalez JR, Carrasco JL, Dudbridge F, Armengol L, Estivill X, Moreno V: Maximizing association statistics over genetic models. Genet Epidemiol. 2008, 32 (3): 246254. 10.1002/gepi.20299.
 15.
Hothorn LA, Hothorn T: Orderrestricted scores test for the evaluation of populationbased case–control studies when the genetic model is unknown. Biometrical J. 2009, 51 (4): 659669. 10.1002/bimj.200800203.
 16.
Purcell S, Neale B, ToddBrown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, de Bakker PIW, Daly MJ, Sham PC: PLINK: a tool set for wholegenome association and populationbased linkage analyses. Am J Hum Genet. 2007, 81 (3): 559575. 10.1086/519795.
 17.
Lu JC, Wei QY, Bondy ML, Yu TK, Li DH, Brewster A, Shete S, Sahin A, MericBernstam F, Wang LE: Promoter polymorphism (−786T > C) in the endothelial nitric oxide synthase gene is associated with risk of sporadic breast cancer in nonHispanic white women age younger than 55 years. Cancer. 2006, 107 (9): 22452253. 10.1002/cncr.22269.
 18.
Casella G, Berger RL: Statistical Inference. 2002, Australia ; Pacific Grove, CA: Thomson Learning, 4654. 2
 19.
Mordell LJ: Diophantine Equations. 1969, Academic P: London, New York
Acknowledgements
This work was supported by National Institutes of Health grants R01CA131324 (SS), NIH R25 DA026120 (SS), and R01DE022891 (SS). This research was supported in part by Barnhart Family Distinguished Professorship in Targeted Therapy (SS). This research was supported in part by a cancer prevention fellowship for Rajesh Talluri supported by a grant from the National Institute of Drug Abuse (NIH R25 DA026120).
Author information
Additional information
Competing interests
We declare that there are no competing interests.
Authors’ contributions
RT and SS conceived and designed the study. RT implemented the method. RT and JW performed simulations. RT and SS wrote the paper. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Genetic association
 Multiple testing
 CochranArmitage trend test
 Genetic models
 Exact pvalue