PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Stat Interface. Author manuscript; available in PMC 2010 December 15.
Published in final edited form as:
Stat Interface. 2009 January 1; 2(2): 161–170.
PMCID: PMC3002050
NIHMSID: NIHMS242896

Detecting essential and removable interactions in genome-wide association studies

Abstract

Detection of disease gene interaction effects among the enormous array of single nucleotide polymorphism (SNP) combinations represents the next frontier in genome-wide association (GWA) studies. Here we propose a novel strategy on the basis of the pattern and nature of the interaction, which can be classified as essential (EI) or removable (RI). We provide an analytical framework, including the qualitative conditions for screening EIs/RIs and a RI-to-EI likelihood ratio score to quantitatively measure the effect. In analyzing six GWA data sets, we find that the scores follow an exponential distribution, except in the upper 10−8 tail region in which the scores become irregular and unpredictable. Our approach is conceptually simple, computationally efficient and detects interactions that can be visualized and unequivocally interpreted.

Keywords: Genome-wide association study, gene-gene interaction, removable interaction, essential interaction, likelihood ratio statistics

1. INTRODUCTION

A growing number of successful genome-wide association (GWA) studies have been reported in the past several years [6]. All have initially focused on the main effects of individual single nucleotide polymorphisms (SNPs), which for the most part have been found to be relatively weak. Yet to be explored in this wealth of data are the interaction effects among SNPs, wherein the main effect of one SNP may be weak because that SNP confers risk only in the presence of one or more other SNPs. Taking these gene-gene interactions into account may lead to the discovery of stronger and more interesting effects [2, 8]. However, searching for meaningful interactions is a challenge for even the simplest case of examining all two-way interactions among a million SNPs (~half of 1012 combinations).

To aid in efficiency and interpretability of this formidable computational problem, we present an analytical framework for two types of interactions, removable interactions (RIs) and essential interactions [1, 12]. Interaction is usually defined as departure from additivity of effects on a specific outcome scale. If a monotone transformation exists that induces additivity, the interaction is called removable; otherwise, the interaction is essential.

What makes the mathematics hard is the sheer volume of networks in GWA data, for which no universal law underlying the system can be deduced. Our idea to detect interactions in GWA studies is based on distinct structures which can be visualized and epidemiological interpretations of the effects. We developed a qualitative tool that can rapidly screen specific interactions, essential or removable, and a likelihood score that quantifies the extent of the effect. We have identified empirical regularities in applying the proposed method to six GWA data sets with number of SNPs ranged from 100,000 to 1,000,000.

Under the simplest scenario of two SNPs (SNP1 and SNP2), each having two genotypic variants denoted by 0 and 1, respectively, we derived the necessary and sufficient conditions to use as a screening tool to quickly and reliably identify pair interactions as essential or removable. We then derived the EI-to-RI likelihood ratio (EI-RI score) to quantify the effects. In investigation of six existing GWA studies we find a consistent pattern of the EI-RI score distributions. Extension to scenarios with SNPs having three genotypes is straightforward, albeit involving tedious calculations.

2. METHODS

2.1 Essential/removable interactions

The joint risks of a pair of dichotomous markers can be graphically represented by two lines and are depicted in Fig. 1. We use the OR and log(OR) risk scales; however the results presented can be generalized to any risk scale. Figure 1 illustrates that all two-SNP combinations can be classified into one of three mutually-exclusive categories:

Figure 1
Graphical representations of absolute non-interaction, removable interaction, and essential interaction. (a) Absolute non-interaction (ANI): in the presence of a given genotypic variant of SNP2, risk does not vary by genotypic variant of SNP1 (two lines ...

Absolute non-interaction (ANI) (Fig. 1a): No interaction between SNP1 and SNP2 means that the effect of SNP1+SNP2 is the sum of the individual effects (i.e., effects are additive). ANI is manifest by a two-SNP combination that does not produce an interaction under any risk scale. This phenomenon occurs when in the presence of a given genotypic variant of SNP2, risk does not vary by genotypic variant of SNP1.

Removable interaction (RI), Fig. 1b–c, is manifest by a two-SNP combination that produces an interaction under at least one risk scale (Fig. 1b), and does not produce an interaction under at least one other risk scale (Fig. 1c). If an RI exists, the direction, positive or negative, of the risk difference for one SNP is not affected by the other SNP (and is not zero). However, under a risk scale in which the RI manifests itself as an interaction, the magnitude of the effect of SNP1 varies by SNP2 genotypic stratum, and vice versa. That is, with an RI, the magnitude, but not the direction, of the effect for one SNP may be modified by the other SNP. An example of an RI is given in Table 1.

Table 1
An example of a removable interaction (RI). The interaction between two SNPs (SNP1 and SNP2) is significant under the odds ratio (OR) scale, but is removed under the log(OR) scale. This example is graphically illustrated in Fig. 1b–c

Essential interaction (EI), Fig. 1d-f, is manifest by a two SNP combination that produces an interaction under all risk scales: In contrast to RI, an EI exists when the direction of the effect of one SNP is dependent on the genotypic variant of the other SNP. That is, the effect of one SNP is reversed by the other SNP. And as in RI, for an EI the magnitude of the effect for one SNP may also be modified by the other SNP. EI can be further classified into two subclasses, EISLOPE and EICROSS. EISLOPE occurs when the reversal of the direction of effect is one way (e.g., SNP1 by SNP2 but not SNP2 by SNP1) (Fig. 1d–e), while EICROSS occurs when the reversal of the direction of effect is reciprocal (SNP1 by SNP2 and SNP2 by SNP1) (Fig. 1f).

From the above reasoning we can derive a set of relationships in terms of magnitudes and directions of the effect change of one SNP by the other:

Condition I (slopes have opposite direction): (OR10OR00)(OR11OR01) ≤ 0 in which one risk difference in the product can equal zero, but not both.

Condition II (lines intersect): (OR01OR00)(OR11OR10) ≤ 0 in which one risk difference, but not both, can equal zero.

An interaction is essential if and only if Condition I and/or II holds. If only one of these conditions holds, we have an EISLOPE; if both conditions hold, we have an EICROSS. If both conditions are not satisfied, then we have RI (unless both product terms in Condition I or Condition II are zeros, in which case we have ANI).

2.2 An EI-RI score: a quantitative measure for the effect

Since the Conditions I-II are invariant under any monotone transformation, ORab in two conditions could be replaced by its logarithm, the log odds ratio of genotypic variants combination ab with respect to baseline 00. Based on the Condition I and/or Condition II, the necessary and sufficient(N&S) condition for EIs, a statistical test can be constructed for the null hypothesis that no EI exists. By replacing ORs with their logarithms, the testing hypothesis depends only on log-odds ratios. The standard theory [11, 13, 14] assures that the likelihood ratio test for such kind of hypothesis based on the prospective likelihood function would be valid. The prospective likelihood function is proportional to

L(θ)=Πa=01Πb=01pabnab(1pab)mab,
(1)

where nab and mab respectively are observed numbers of cases and controls for joint genotypic variant ab, pab = eθab /(1 + eθab) is the corresponding affected probability and θab represents the log odd of the joint genotypic variant ab. Let θ = (θ00, θ01, θ10, θ11) be the vector of parameters. The null hypothesis parameter space, the complement of the subset for parameters satisfying the N&S condition, is not closed. We will use instead its closure, i.e.,

ϴ0={θ:(θ01θ00)(θ11θ10)0and(θ01θ00)(θ11θ10)0},
(2)

as a working null parameter space. By doing so, we can simplify the computation and the test will be slightly conservative.

We took the EI-RI score to be the likelihood ratio statistics [7], i.e.

2(maxθϴlogL(θ)maxθϴlogL(θ)),
(3)

where ϴ = {θ : −∞ ≤ θ00, θ01, θ10, θ11 < is the uncon-strained parameter space.

Obviously, the unconstrained MLE (θ [set membership] ϴ) has a closed form, namely θab = log(nab/mab). To obtain the EI-RI score, one needs to calculate the constrained MLE under the null hypothesis, which can be solved by a non-linear programming (NLP) algorithm. However, the NLP algorithm is time-consuming and could break down especially when the constrained MLE is on the boundary of ϴ0. Therefore, the NLP algorithm is not suitable for GWA study.

By virtue of the concavity of the log likelihood function, we developed an efficient algorithm for calculating the constrained MLE. Details of the algorithm are given in Appendix A.

2.3 EI-RI scores in six GWA data sets

Under the null hypothesis of ANI or RI, the null distribution of the EI-RI score is mathematically intractable as it is a mixture of zero (when the maximum point belongs to the null parameter space) and a positive distribution (when the constrained MLE does not belong to the null parameter space) with the component weights unknown. Estimation of p-values by the permutation test fails because the conventional permutation procedure is not suitable for our composite null hypothesis. In fact, randomly allocating the samples under the assumption of no joint effects or marginal effects as in the conventional permutation approach inevitably introduces false positives. And indeed, serious inflated type I error rates by the permutation test have been observed in our simulation studies (data not shown here).

Re-sampling approaches [4] might be applicable for estimating p-values with small numbers of SNP markers, but it is beyond computational capability for GWA data with large numbers of SNP pairs. We therefore sought to examine the empirical distribution patterns of the scores among six GWA data sets. All empirical distributions of EI-RI scores in six data sets give rise to similar upper tails, which, in turn, may justify the existence of a universal threshold for declaring the significance of an EI.

3. RESULTS FOR EMPIRICAL STUDY

Descriptive statistics of six GWA data sets are given in Table 2. The numbers of SNPs in these studies ranged from 105 to 106. We dichotomized each SNP genotype based on the empirical mode of inheritance [9]. In the six data sets examined, about half (48.83% on average) of the SNP pairs were in the category of EIs (Fig. 1d–f). Of the EI pairs, about two-thirds (65.97% on average) were EISLOPE pairs and one-third (34.03%) were EICROSS pairs.

Table 2
Descriptive information of six GWA data sets

Histogram plots are drawn for EI-RI scores of the six data sets and are shown in Fig. 2a. Comparing histogram plots of the two EI groups, EICROSS pairs tend to have higher EI-RI scores than EISLOPE pairs. Patterns of histogram plots are consistent across six data sets, particularly, when EI-RI scores become large. This consistency implies that the empirical distributions of EI pairs in the six data sets have similar tail behavior. This conclusion was confirmed by the consistency of tail quantiles for EI-RI scores, as shown in Fig. 2c. The consistent tail distribution of EI-RI scores in six data sets may justify the existence of a universal critical value for declaring the significance of an EI.

Figure 2a
The histogram plots of EI-RI scores in six data sets. Dark black is the histogram for EICROSS pairs and gray is the histogram for EISLOPE pairs. For each type of EI pair, the leftmost column is the frequency for EI-RI scores less than 0.0005 and the rightmost ...
Figure 2c
The tail quantiles of the EI-RI scores stratified by EICROSS pairs and EISLOPE pairs in the six GWA data sets based on the real data (solid lines) and the extrapolated tail quantiles from an estimated exponential distribution (dotted lines).

In GWA studies, the critical values are chosen as the (1-p)-quantiles of EI-RI scores under the null hypothesis for an extremely small p. One way to estimate the extreme quantiles is to inverse the linear interpolation of the empirical cumulative distribution function [10]. We found that the empirical (1-p)-quantiles are consistent across the six data sets for 10−8p ≤ 10−2 (Fig. 2c, Table 3). However, in the upper 10−8 tail region the quantiles became irregular and unpredictable, which is due to random error of the extreme value distribution of p-values. Therefore, the interpolation method is not suitable for estimation of extreme quantiles. To overcome this problem, we approximated extreme quantiles using much lower quantiles which can be estimated by intermediate order statistics.

Table 3
Right tail quantiles of the EI scores in the six GWA data sets. Quantiles between 10−8 and 10−2 were estimated directly from the data, outside that range were estimated by extrapolation

For each empirical distribution of the six data sets, the relationship between the (1 – p) quantile and log(p) was approximately linear in the range 10−8p ≤ 10−2 (Fig. 2c). The linearity implies that extreme value indices [5] are zero for all six empirical distributions. The conclusion was confirmed by constructing 95% confidence intervals of the extreme value index with different numbers of upper orderstatistics (for details see Appendix B). The value zero belongs to almost all these confidence intervals, which does not contradict the hypothesis that the extreme value index is zero.

With a zero extreme value index, the (1 − p)-quantile for small p can be approximated by

Q(1p)=σlogp+μ
(4)

for some μ and σ(σ > 0). That is, the upper tail of the EI-RI score approximately follows an exponential distribution with location parameter σ and scale parameter σ.The two parameters were estimated by fitting the regression line defined in (4) for p in the range [10−8, 10−2]. As shown in Table 4, all R-square statistics are close to 1, which indicates that regression lines defined in (4) almost perfectly fit the data (also see quantile-quantile plots in Fig. 2b).

Figure 2b
Quantile-quantile plots of EI-RI scores for six GWA data sets. Each point corresponds to a probability p: the x-coordinate represents the EI-RI score corresponding to the p-th quantile of the exponential distribution (details in text) and the y-coordinate ...
Table 4
Estimation for intercept μ and slope σ by fitting (1 – p)-empirical quantile vs −log(p) for 10−8 ≤ p ≤ 10−2

Furthermore, estimators of the two parameters are consistent across all six data sets. This, again, suggests that distributions of EI-RI scores have similar tail behavior.

Based on these data and a reasonable underlying assumption of an exponential tail distribution for EI-RI scores under the null hypothesis, we can deduce a simple proposition: the EI-RI score corresponding to the (1 − p)-quantile is approximately equal to (−2.617) × log(p) – 2.483, for small p, and equivalently, the p-value is approximately equal to exp{−(EI-RI score + 2.483)/2.617}, for large EI-RI scores.

It will be interesting to see whether a similar trend of the scores will hold and the above estimated p-values can be generalized for any GWA study. If so, this in turn may provide a simple guideline of choosing significance levels for SNP pairs worthy of further investigation such as replication studies. For example, here we postulate an EI-RI score less than 18 indicates an inconclusive interaction; a score between 18 and 26 suggests a possibility of a true EI; and a score greater than 26 indicates a statistically significant EI. These speculations warrant further investigation.

4. DISCUSSION

Detection of interaction effects among SNPs represents the next frontier in GWA studies. However, screening the enormous array of SNP combinations represents a daunting task. To begin to address this challenge, we classified all two-SNP combinations into three mutually-exclusive categories, absolute non-interactions, removable interactions, and essential interactions, and we characterized the conditions for each of these categories. In essence, an interaction is removable when the magnitude, but not the direction, of the effect of one SNP is modified by the other SNP, and an interaction is essential when the direction of the effect for at least one of the SNPs is altered in the presence of the other SNP (and the magnitude may be modified as well).

We developed an efficient computer program to screen for EIs and to calculate the EI-RI scores (likelihood ratio statistics) for each SNP pair in GWA data of any number of SNPs. With this algorithm we examined the distributions of EI-RI scores in six GWA data sets with 105 to 106 SNPs. A consistent pattern emerged, which allowed us to derive the tail quantiles by extrapolation from the stable portion of the empirical distributions. Finally, we suggested an EI-RI score threshold for identifying prime candidate EIs.

The EI-RI score presented here emphasizes EIs because they are not dependent on risk scale and represent the most extreme form of interaction (reversal of direction of risk difference). Furthermore, while one may reasonably screen for RIs solely among SNPs that exhibit significant main effects, this approach may not work for EIs, especially for EICROSS. This is because in the presence of an RI, both SNPs have non-zero marginal effects (of course, the magnitude of the effects depends on the joint effects and MAFs). With EI, however, the marginal effect for at least one SNP can be close to zero even when the EI effect is substantial. With EICROSS the marginal effect for both SNPs can be close to zero. Thus, many EIs may be missed if only SNPs with main effects are screened. This conjecture requires further investigation.

Even though RIs have the undesirable mathematical property of dependence on the risk scale, it is likely that many RIs in GWA data sets will prove to be potentially interesting, as they have proved to be in traditional case-control studies. It will be straightforward to apply the methods presented here, which are independent of risk scale, to screen for RIs as well as EIs.

Taken together, the qualitative and the quantitative methods we proposed to detect interaction effects in GWA data by way of EIs and RIs is conceptually simple, computationally manageable and, most of all, the resulting interactions can be readily visualized and unequivocally interpreted. We plan to determine the practical utility of our current methodology to address real disease problems.

All algorithms have been implemented in C++ and the program is available upon request from authors.

ACKNOWLEDGEMENTS

The research was partially supported by NIH/NEI grants (R01 EY015771, R21 EY018127), NIH/NHGRI grant (K25 HG000060), NIH grant (5-R37-GM047845-18), NSF grant (DMS-05-04871), NSFC grants (10701067, 10671189), and grants provided by Verto Institute, American Health Alliance Foundation, Macular Vision Research Foundation, and Ellison Medical Researching Foundation. This work was carried out while Zhang was a postdoctoral associate and Yang was a visiting scholar in Hoh’s Lab at Yale. We gratefully acknowledge Yale University Biomedical High Performance Computing Center and NIH grant: RR19895, which funded the instrumentation. Our thanks also extend to the database of Genotypes and Phenotypes (dbGaP) for the GWA data sets. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors would like to thank the editor for his comments which improved this article.

APPENDIX A. ALGORITHM FOR CONSTRAINED MLE

If the unconstrained MLE belongs to ϴ0, then it is of course the constrained MLE. If the constrained MLE is an interior point of ϴ0, it must be the unconstrained MLE since the log likelihood function is concave so that there are no local maxima other than the global maxima. Therefore, if the unconstrained MLE does not belong to ϴ0, the constrained MLE must be on its boundary ϴB. This leads to the algorithm for obtaining the constrained MLE, which is illustrated in Fig. A1.

Figure A
Algorithm for calculating MLE constrained on the null hypothesis. (1) Calculation of constrained MLE θ^c. (2) Calculation of constrained MLE θ^B1 on boundary ϴB1.

The boundary ϴB is the union of four subboundaries: ϴB1={θ:θ01=θ00and(θ10θ00)(θ11θ01)0}, ϴB2={θ:θ10=θ11and(θ10θ00)(θ11θ01)0}, ϴB3={θ:θ10=θ00and(θ01θ00)(θ11θ10)0} and ϴB4={θ:θ11=θ01and(θ01θ00)(θ11θ01)0}. We can calculate the MLEs on four sub-boundaries separately and the one with largest likelihood function value is the constrained MLE on boundary ϴB. We also develop an algorithm for obtaining the MLE on sub-boundary ϴB1(the other three constrained MLEs can be similarly obtained). The algorithm is graphically displayed in Fig. A2.

The algorithm is very easy to implement and is time efficient since there are closed forms for these constrained MLEs that need to be calculated in the algorithm. For example, under constraint θ01 = θ00, the resulting MLE of θ is (log[(n00 + n01)/(m00 + m01)], log[(n00 + n01)/(m00 + m01)], log[n10/m10], log[n11/m11]). Under constraint θ00 = θ01 = θ10, MLEs of θijs are θ00 = θ01 = θ10 = log[(n00 + n01 + n10)/(m00 + m01 + m10)] and θ11 = log[n11/m11]. Similarly, under constraint θ00 = θ01 = θ11, the MLE of θ00 = θ01 = θ11 is log[(n00 + n01 + n11)/(m00 + m01 + m11)] and the MLE of θ10 is log[n10/m10].

APPENDIX B. MOMENT ESTIMATION AND CONFIDENCE INTERVAL FOR EXTREME VALUE INDEX

The key issue in extreme value theory is the estimation of the extreme value index c, which governs the tail behavior of a distribution function. The Hill estimator is the most popular estimator which is restricted to the case c > 0. The moment estimator c^M, an adaptation of the Hill estimator, is a consistent estimator for all real value c’s. Let X1,nX2, n(...)Xn,n be the order statistics for sample X1, … , Xn, the moment estimator of c is given by

c^M=M1+112(1M12M2)1,

with Mj=1kΣk=0k1(logXni,nlogXnk,n)j for j = 1,2.

Under some regular conditions, k(c^Mc) asymptotically follows a zero mean normal distribution with variance [3]

{1+c2,ifc0(1c)2(12c)(1c+6c2)(13c)(14c),ifc<0.}

An approximate (1 – α)100% confidence interval is then given by

c^MZα2var(c^M)k<c<c^M+Zα2var(c^M)k

where var(c^M) is the asymptotic variance with c being replaced by its moment estimator, and Zα/2 is the (1 − α/2)-quantile of the standard normal distribution.

Table A gives the 95% asymptotic confidence intervals for various values of k. The value zero belongs to almost all these confidence intervals, which does not contradict the hypothesis that the extreme value index is zero.

Table A
Moment estimator and the 95% asymptotic confidence interval for extreme value index c with different k’s, where k is the number of upper order statistics used in the moment estimator

Contributor Information

Chengqing Wu, Yale School of Public Health, New Haven, CT, ude.elay@uw.gniqgnehc.

Hong Zhang, Department of Statistics and Finance, University of Science and Technology of China, Hefei, Anhui, P. R. China, nc.ude.ctsu@hgnahz.

Xiangtao Liu, Department of Applied Mathematics, Yale University, New Haven, CT, ude.elay@uil.oatgnaix.

Andrew DeWan, Yale School of Public Health, New Haven, CT, ude.elay@nawed.werdna.

Robert Dubrow, Yale School of Public Health, New Haven, CT, ude.elay@worbud.trebor.

Zhiliang Ying, Department of Statistics, Columbia University, New York, NY, ude.aibmuloc.tats@gniyz.

Yaning Yang, Department of Statistics and Finance, University of Science and Technology of China, Hefei, Anhui, P. R. China, nc.ude.ctsu@gnayny.

Josephine Hoh, Yale School of Public Health, New Haven, CT, ude.elay@hoh.enihpesoj.

REFERENCES

[1] Breslow NE, Day NE. Statistical Methods in Cancer Research. Volume I - The Analysis of Case-Control Studies. IARC Sci. Publ.; Lyon, France: 1980. [PubMed]
[2] Chen X, Liu CT, Zhang M, Zhang H. A forest-based approach to identifying gene and gene gene interactions. Proc Natl Acad Sci. 2007;104:19199–19203. [PubMed]
[3] Dekkers ALM, Einmahl JHJ, Haan LD. A moment estimator for the index of an extreme-value distribution. The Annals of Statistics. 1989;17:1833–1855. MR1026315.
[4] Efron B. The Jackknife, the Bootstrap and Other Resampling Plans. SIAM; Philadelphia: 1982. MR0659849.
[5] Haan LD, Ferreira A. Extreme Value Theory: An Introduction. Springer; New York, London: 2006. MR2234156.
[6] Hindorff L, Junkins H, Manolio T. A Catalog of Published Genome-Wide Association Studies. 2009. Available at: www.genome.gov/26525384.
[7] Lehmann EL, Romano JP. Testing Statistical Hypotheses. Springer; New York: 2005. MR2135927.
[8] Marchini J, Donnelly P, Cardon LR. Genome-wide strategies for detecting multiple loci that influence complex diseases. Nat Genet. 2005;37:413–417. [PubMed]
[9] Matthews AG, Haynes C, Liu C, Ott J. Collapsing SNP genotypes in case-control genome-wide association studies increases the type I error rate and power. Stat Appl Genet Mol Biol. 2008;7 Article23. [PMC free article] [PubMed]
[10] Parzen E. Nonparametric statistical data modeling. Journal of the American Statistical Association. 1979;74:105–121. MR0529528.
[11] Prentice RL, Pyke R. Logistic disease incidence models and case-control studies. Biometrika. 1979;66:403–411. MR0556730.
[12] Scheffe H. The Analysis of Variance. John Wiley & Sons; New York: 1959. MR0116429.
[13] Scott AJ, Wild CJ. Hypothesis testing in case-control studies. Biometrika. 1989;76:806–808. MR1041428.
[14] Seaman SR, Richardson S. Equivalence of prospective and retrospective models in the Bayesian analysis of case-control studies. Biometrika. 2004;91:15–25. MR2050457.