PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of hheKargerHomeAlertsResources
 
Hum Hered. 2009 October; 69(1): 14–27.
Published online 2009 October 2. doi:  10.1159/000243150
PMCID: PMC2880732

Allelic Based Gene-Gene Interaction in Case-Control Studies

Abstract

In case-control studies identifying disease susceptibility loci, it has been shown that the interaction caused by multiple single nucleotide polymorphisms (SNPs) within a gene as well as by SNPs at unlinked genes plays an important role in influencing risk of a disease. A novel statistical approach is proposed to detect gene-gene interactions at the allelic level contributing to a disease trait. With a new allelic score inferred from the observed genotypes at two or more unlinked SNPs, we derive a score test from logistic regression and test for association of the allelic scores with a disease trait. Furthermore, F and likelihood ratio tests are derived from Cochran-Armitage regression. By testing for the association, the interaction can be assessed both in cases where the SNP association can be detected and cannot be detected as a main effect in single SNP approach. The analytical power and type I error rates over 6 two-way interaction models are investigated based on the non-centrality parameter approximation of the score test. Simulation studies demonstrate that (1) the power of the score test is asymptotically equivalent to that of the test statistics by the Cochran-Armitage method and (2) the allelic based method provides higher power than two genotypic based methods.

Key Words: Allelic test, Interaction effect, Score test, Cochran-Armitage method, Epistasis

Introduction

With the advance of high-throughput sequencing technologies, genetic association studies to identify susceptibility genes of a common, complex, human disease has been promoted. It has been known that multiple genes along with the environment play interactive roles in contribution to the development of complex diseases. Most ongoing genetic studies have focused on identifying the effect of a single gene; these studies often fail to identify causal association with a disease due to the fact that a disease is the consequence of a complicated network of multiple susceptibility loci, each of which is likely to have only a small effect when considered alone.

Recently many statistical methods have been proposed to detect gene-gene interactions in a case-control study. As a machine learning and data mining approach, multifactor-dimensionality reduction (MDR) was proposed to pool multilocus genotypes into high-risk and low-risk groups and analyze the datasets through cross-validation to calculate prediction accuracy of the testing SNPs [1]. Classical decision tree approaches such as classification and regression trees [2] and support vector machines [3] classify subjects into case and control groups by an optimal hyperplane predicted by SNPs combination. Similar to other statistical model-based methods, the logit regression method [4] is a generalized regression model that uses a ranking criterion to identify the best model. The most commonly used method is a multiplicative method in a standard logistic model framework which consists of main effects and multiplicative interaction effects of SNPs.

Most currently available methods of either machine learning approaches or model based approaches focus on the association of genotypic combinations of multiple SNPs with a disease trait. Jung et al. [5] proposed a novel approach to search for the interaction at the allelic level for quantitative trait loci (QTL). As stated in Jung et al. [5], there are several differences in comparison between genotype-based methods and allelic-based method. Genotypic-based approaches can increase the degrees of freedom of a test statistic, resulting in a concomitant loss of power when multiple unlinked markers (i.e. more than 3) are tested for interaction. In contrast, the allelic-based method can reduce the number of degrees of freedom by accounting for the allelic levels. Furthermore, the power of the genotypic-based method may be compromised by sparse or empty combination caused by multiple markers with small sample sizes or a lower minor allele frequency (≤ 0.1). However, the allelic-based method is relatively robust even with the lower minor allele frequency and small sample sizes.

In this paper, we extend our previous work detecting interaction of quantitative trait loci (QTL) to a case-control study design comparing those with and those without the disease. By a new definition of the allelic-based gene-gene interaction, which is the nonrandom association of alleles occurring when a particular allele in one gene and a particular allele in another unlinked gene contribute to the risk of a disease trait through their interaction [6], the proposed method tests for association of these allelic scores assigned to each subject with a disease trait. We derive a score test based on the standard logistic regression model and both F test and likelihood ratio test (LRT) statistics based on the Cochran-Armitage (CA) regression model. In addition, the non-centrality parameters approximation of these statistics were derived to evaluate the analytical properties of the test statistics. We perform simulation studies to demonstrate the analytical properties of the score test and CA-based tests and to compare the power of the allelic-based method with that of the genotypic-based methods. The 7 candidate genes of major depression disorder (MDD) case-control data provided in dbGaP [7] were applied to show the feasibilities of the proposed method.

Methods

The allelic-based gene-gene approach proposes to test for non-random association of the allelic combination with a disease trait. The allelic combination of each subject is scored by a conditional probability of having the particular allelic combination given the observed genotypes at the putative unlinked loci of each subject.

Statistical Models and Properties

Let D*1, D*2 be two disease loci underlying a disease trait, which are both in Hardy-Weinberg equilibrium (HWE) and are unlinked but associated with a disease through their interaction. Let D1 and d1 be two alleles at the first disease locus D*1, with frequencies pD1 and pd1 respectively. Let D2 and d2 be two alleles at the second disease locus D*2, with frequencies pD2 and pd2. We consider two observed marker loci, each of which is in linkage disequilibrium (LD) with either of the two interacting disease loci; Assume that marker M1 is in LD with the susceptibility locus D*1 and marker M2 is in LD with the other susceptibility locus D*2 and both are unlinked. An observed marker M1 has two alleles A and a with frequencies pA and pa respectively, and M2 has two alleles B and b with frequencies pB and pb respectively. Let DAD1 = P (AD1) – PAPD1 be the measure of LD between D*1 and M1, and DBD2 = P (BD2) – PBPD2 be LD measure between D*2 and M2. Denote the penetrance of disease genotype of D*1 by fkl = p (Y = 1 | D*1 [sm epsilon] (D1D1, D1d1, d1d1)), k, l [sm epsilon] (D1, d1), and that of D*2 by fmn = p (Y = 1 | D*2 [sm epsilon] (D2D2, D2d2, d2d2)), m, n [sm epsilon] (D2, d2) respectively. Let fklmn = P (Y = 1 | D*1, D*2) be the joint penetrance for the two disease loci D*1 and D*2 k, l [sm epsilon] (D1, d1) m, n [sm epsilon] (D2, d2). As shown in figure figure1,1, the penetrance of a particular allelic combination is calculated using the joint penetrance of genotypic combinations consisting of the particular alleles. Therefore, we can derive the penetrance of two D1D2 allelic combination by fD1D2 = P (Y = 1 | D1D2) = fD1D1D2D2 pD1 pD2 + fD1D1D2d2 pD1 pd2 + fD1d1D2D2 pd1 pD2 + fD1d1D2d2 pd1 pd2. In the same way fD1d2, fd1D2, fd1d2 can also be calculated. For a control group, the penetrance f is replaced by 1 – f such as 1 – fkl = p (Y = 0 | D*1 [sm epsilon] (D1, D1, D1d1, d1d1)) and 1 – fmn = p (Y = 0 | D*2 [sm epsilon] (D2, D2, D2d2, d2d2)), k, l [sm epsilon] (D1, d1), m, n [sm epsilon] (D2, d2) for the two genes.

Fig. 1
Cell combinations of unlinked two disease loci.

Logistic Regression Model and Non-Centrality Parameter

Denote yi = 1 if the i-th subject is affected by a disease (case) and yi = 0 if the i-th subject is not affected by the disease (control). In the non-parametric maximum likelihood solution that allows an arbitrary covariate distribution, Anderson [8], Prentice and Pyke [9] and Scott and Wild [10] showed that fitting a standard prospective logistic regression in a case-control sampling design is equivalent to fitting a retrospective logistic regression except that an intercept in the case-control sampling needs the information of sampling fraction of cases and controls. Therefore, the standard logistic regression model is used due to the equivalence in the parameter estimates of interaction effect. The likelihood function of the standard logistic regression is

L(α,βAB,βAb,βaB)=i=1N{[P(yi=1|Xi,AB,Xi,aB,Xi,Ab,wi)]yi×[1-P(yi=1|Xi,AB,Xi,aB,Xi,Ab,wi)]1-yi},
(1)

where

P(yi=1|Xi,AB,Xi,aB,Xi,Ab,wi)=exp(α+w1γ+Σk={AB,aB,Ab}Xi,kβk)/[1+exp(α+wiγ+Σk={AB,aB,Ab}Xi,kβk)].

Xi,AB = P(AB | M1, M2) is an allelic score of AB allelic combination from M1 and M2 genotypes of the i-th subject and βk [sm epsilon] {AB,Ab,aB} is the interaction effects of the k-th [sm epsilon] { AB, Ab, aB } allelic combination. wi is a covariate such as age or gender and is a coefficient of the covariate. The overall proportion of y is R / N, where R is the number of case subjects and N is the number of total subjects. As in Jung et al. [5], table table11 presents the allelic scores given the observed genotypes of two markers. The number of allelic combinations at two markers is 22 but since the probability of a particular allele combination is a function of the rest of the probability of three probabilities, so the number of independent variables is 22 − 1. The null hypothesis (H0 : βAB = βAb = βaB = 0) is that there is no association caused by the interaction; therefore, disease risk of the four allelic combinations are not significantly different. Rejecting the null hypothesis would indicate there is an allelic interaction between two markers, especially particular allelic combinations generated by the different genetic characteristics between cases and controls.

Table 1
Allelic scores given genotypes of individual

Under the assumption of no covariates from (1), an intercept is a nuisance parameter. Let Uτ = (UAB, UAb, UaB)τ be a score function, which is the derivative of the log-likelihood function with respect to β = (βAB, βAb, βaB) respectively. Under the null hypothesis H0 : βAB = βAb = βaB = 0 of no interaction effect, the maximum likelihood estimate (MLE) of the intercept is αˆ=log(y¯/1-y¯). The efficient score test statistic under the null hypothesis is

S=UτV-1U,
(2)

where

Uτ=Uβ(βH0,αˆ)=(Σxi,AB(yi-y¯),Σxi,Ab(yi-y¯),Σxi,aB(yi-y¯))τ

and V− 1 is the submatrix of I–1 (α, βAB, βAb, βaB) corresponding to Uτ = (UAB, UAb, UaB)τ. Appendix A shows the detailed derivation of Uτ and I (α, βAB, βAb, βaB). The test statistic S is a central chi-square distribution on 3 degrees of freedom.

To investigate the analytical power of the score test, we need to approximate a distribution of S under the alternative hypothesis [11]. The approximation to U under the alternative hypothesis is

U*=Uβ(βH0,α*)-Iβα*Iαα*-1Uα(βH0,α*)=Σi(yi-y¯)[(Xi,ABXi,AbXi,aB)-(x¯ABx¯Abx¯aB)],
(3)

where UαH0, α*) = [partial differential]log L/[partial differential]α, I*βα, I*αα are a submatrix of the Fisher information matrix, I, respectively (see Appendix A) and α* is a solution of the equation 0 = lim n−1E [UαH0,α)]. Let Ω and Σ denote the mean and covariance matrix of U* [11]. The distribution of S* = U Σ−1U* is approximately a chi-square distribution on 3 degrees of freedom with non-centrality parameter,

λscore=ΩτΣ-1Ω,

where

Ω=E(U*)=[(fD1D2-fD1d2-fd1D2-fd1d2)(DAD1DBD2DAD1DbD2DaD1DBD2)+(fD1-fd1)(DAD1PBDAD1PbDaD1PB)+(fD2-fd2)(DBD2PADbD2PADBD2Pa)],

Σ=V(U*)=NΦ(1-Φ)[(PAPB(1-0.5Pa)(1-0.5Pb)PAPBPb(1-0.5Pa)PAPBPa(1-0.5Pb)PAPBPb(1-0.5Pa)PAPb(1-0.5Pa)(1-0.5PB)PAPBPaPbPAPBPa(1-0.5Pb)PAPBPaPbPaPB(1-0.5PA)(1-0.5Pb))-Φ(1-Φ)PXtPX],

where Φ = R / N. Furthermore, (fD1fd1), (fD2fd2) can be explained by the average effect of the gene substitution at each D*1, D*2, [11]. (fD1D2fD1d2fd1D2 + fd1d2) is the magnitude of interaction effect. The non-centrality parameter is also a function of LD between a disease locus and a marker. Appendix B shows the detailed derivation of Ω and Σ. We can further investigate which allelic combinations are interacting to influence a disease trait by testing the allelic specific interaction effects under the null hypothesis of H0 : βC = 0, in the logistic model with the allelic combinations of interest, logit[P(yi = 1 | wi, Xi,C)] = α + wiγ + Xi,C βC. In the extension to multiple markers, the assignment of allelic scores to each subject and modeling logistic regression are straightforward.

Cochran-Armitage Trend Regression Model and Non-Centrality Parameter

With the same allelic scores at two markers in table table1,1, we can model a linear trend of proportions of cases to total (sum of case and control) at each allelic combination, pk,j = rk,j / nk, j, where nk,j = rk,j + sk, j for k [sm epsilon] (AB, Ab, aB, ab), j [sm epsilon] (0,1/4,1/2,1). rk,j, sk,j are the number of affected subjects and unaffected subjects having j score at k allelic combination respectively. Table Table22 shows the distribution of the allelic combination AB, ZAB,j with the order of (ZAB,0 < ZAB,1/4 < … < ZAB,1) in case-control data and the remaining combinations ZAb,j, ZaB,j can be expressed in the same way. It has been shown that regressing pk,j on ZAB,j, ZAb,j, ZaB,j is equivalent to regressing yi on ZAB,j, ZAb,j, ZaB,j [13, 14].

Table 2
Distribution of allelic combination, AB for case-control data

As an extension of CA method, the interaction effect of two markers can be modeled as

yi=α+wiγ+ZAB,iβAB+ZAb,iβAb+ZaB,iβaB+ε,
(4)

where wi is a covariate such as age or gender and γ is the coefficient of covariate. Under the assumption of no covariates for simplicity, the theoretical regression coefficients can be derived by functions of the genetic effects and of linkage disequilibrium (LD) between a marker and a disease susceptibility locus as follows.

β=(αβABβAbβaB)=K-1[RN(1PAPBPAPbPaPB)+(fD1D2-fD1d2-fd1D2+fd1d2)(0DAD1DBD2DAD1DbD2DaD1DBD2)+(fD1-fd1)(0DAD1PBDAD1PbDaD1PB)+(fD2-fd2)(0DBD2PADbD2PADBD2Pa)],

where

K=(1PAPBPAPbPaPBPAPBPAPB(1-0.5Pa)(1-0.5Pb)0.5PAPBPb(1-0.5Pa)0.5PAPBPa(1-0.5Pb)PAPb0.5PAPBPb(1-0.5Pa)PAPb(1-0.5Pa)(1-0.5PB)0.25PAPBPaPbPaPB0.5PAPBPa(1-0.5Pb)0.25PAPBPaPbPaPB(1-0.5PA)(1-0.5Pb)).

Likewise, (fD1fd1), fD2fd2) are the average effect of the gene substitution at each D*1, D*2 and (fD1D2fD1d2fd1D2 + fd1d2) is the magnitude of interaction effect. Appendix C shows the detailed derivation of the theoretical regression coefficients of the model (4). The global test statistic for the interaction effect is

FCA=(Hβˆ)t[H(XtX)-1Ht]-1(Hβˆ)Yt[IN-X(XtX)-1Xt]YN-43,

where IN is the N × N identity matrix and a test matrix is

H=(010000100001).

The FCA follows F (3, N − 4) with non-centrality parameter λCA = 0 under the null hypothesis H0 : βAB = βAb = βaB = 0. Under the alternative hypothesis, FCA is noncentral to F(3, N– 4) with its non-centrality parameter,

λCAR+Sσ2(βAB,βAb,βaB)[HE(XtX)-1Ht]-1(βAB,βAb,βaB)t,

where σ2 is the total residual variance and N = R + S. Furthermore, an LRT can be derived for the global interaction test of H0 : βAB = βaB = βAb = 0, which is 2(LHALH0) following χ2df = 3. LHA is a likelihood function of maximum likelihood estimates (MLE) αˆ,βˆAB,βˆAb,βˆaB under the alternative hypothesis and LH0 is a likelihood function of MLE of αˆ under the null hypothesis. Likewise, a specific allelic interaction effect under the null hypothesis of H0 : βC = 0 for yi = α + wi γ + Xi,CβC can be tested in the same framework.

Model Comparisons

Due to the possibility that the significant allelic interaction results from the additive effect of one or both SNPs, the allelic based model needs to be compared to a main effects model in order to distinguish the real allelic interaction from the main effects. First, we compare the two-way interaction model with a main effect model consisting only of main effects of two putative SNPs. When a logistic regression model is utilized, the criteria for the model comparison are Akaike Information Criterion (AIC) and the smallest Bayesian Information Criterion (BIC) because of the non-nested terms of an interaction model and a main effect model. When CA regression is utilized, an artificial nesting approach called a J test [15] can be employed as follows:

yi=(1-α)f(X,β)+αg(Z,δˆ),

where f(X,β)=μ+Xi,ABβAB+Xi,AbβAb+Xi,aBβaB is a two-way interaction model and g(X,δ) = μ + δAZM1 + δBZM2 is a main effect model. Since f is linear in β, the comparison requires that one estimates and then fits a linear regression and test for α = 0 using the ordinary t-statistic [15]. In addition to the test of α, AIC and BIC are also the criteria to select the best models. The same criteria are applied to the three-way interaction model comparison.

Results

Type I Error Rates

To evaluate the statistical robustness of the score test, the F test and the LRT for the allelic based method, simulation studies were carried out to estimate type I error rates at the 1%, 5% significance level. Six models of interaction between two unlinked disease loci were considered (see table table3).3). Most of the models were designed based on combinations of dominant and recessive inheritance at the genotypic level at each marker. These models are (1) Dominant or Recessive (Dom [dbl union] Rec); (2) Recessive or Recessive (Rec [dbl union] Rec); (3) Dominant and Dominant (Dom ∩ Dom); (4) Dominant and Recessive (Dom ∩ Rec); (5) Threshold whose disease risk is increased when two or more high risk alleles from either locus are present, and (6) Modified model whose homozygosity at either locus confers disease risk. For each model, we simulated 5,000 datasets using SNaP [16]. Each dataset has 200 cases and 200 controls with two unlinked genes under no LD between markers and disease loci (DAD1 = DAD1 / Dmax = 0; DBD2 = DBD2 / Dmax = 0; R2AD1 = R2BD2 = 0), disease risk allele frequencies at two loci of 0.2 (PD1 = PD2 = 0.2), and minor allele frequencies at M1 and M2 of 0.3 (PA = PB = 0.3). Table Table44 shows the results of the empirical type I error rates of the score test and the F and LRT statistic at 1 and 5% significance level. It illustrates that all type I error rates of the three statistics at each model were close to the nominal values 1% and 5% which suggest that the proposed method is statistically robust. For different parameters such as 0.5 disease-associated allele frequencies (PD1 = PD2 = 0.5) and equal allelic frequencies of the two markers (PA = PB = 0.5), all models were also close to nominal values of 1 and 5% as in table table44 (results not shown).

Table 3
The interaction models: f stands for penetrance of genotypic combination
Table 4
Result of type I error rate (%) at 1% and 5% significance level from 5,000 datasets, each has 200 cases and 200 controls

Analytical Power and Sample Size Calculation

The analytical power and the required sample size of the proposed approach were calculated based on the non-centrality parameter λscore = ΩτΣ−1Ω of the score test. The power to detect the interaction effect was influenced by various parameters such as LD between the putative disease loci and a marker (DAD1, DBD2), disease allele frequencies (pD1, pD2), observed marker allele frequencies (pA, pB), and penetrance (f) which is a proportion of individuals carrying a genotypic combination of two putative disease loci that also develop a disease. For the six two-loci interaction models, the power is calculated against each parameter respectively while the other remaining parameters are fixed at a number that provides the greatest variation across the interaction models.

Figure Figure22 shows the power curve of the score test S plotted against penetrance f (see table table3)3) at the 1% significance level. Power is calculated when sample size is 1,000 cases and 1,000 controls with a fixed value DAD1 = DBD2 = 0.05, PA = Pa = 0.5; PB = Pb = 0.5. Models ‘Dom ∩ Dom’ and ‘Threshold’ equivalently achieved the highest power and ‘Modified’ obtained the lowest power across the full range of penetrance f. The power of ‘Rec [dbl union] Rec’ is higher than that of ‘Dom ∩ Rec’ which in turn is higher than that of ‘Dom [dbl union] Rec’ model. Figure Figure33 presents the power against the A allele frequency (pA) of a marker M1, when pB = 0.5 at a marker M2. All parameters used in figure figure33 are the same as those used in figure figure22 except DAD1 = {min(PD1, PA) – PD1PA}/2 and the penetrance f = 0.5 at each model. The pattern of the power in figure figure33 is very similar to that of the power in figure figure22 across the six models except ‘Dom [dbl intersection] Rec’ which achieves the second highest power when pA > 0.2, whereas ‘Rec [dbl union] Rec’ is the second highest power across the full range of penetrance in figure figure22.

Fig. 2
Analytical power against penetrance f at 1% significance level. PD1 = 0.5; PD2 = 0.5; PA = 0.5; PB = 0.5; DAD2 = DBD2 = 0.05; the sample size of cases R = 1,000 and controls S = 1,000. Dominant or Recessive model is indicated by ‘Dom [dbl union] ...
Fig. 3
Analytical power against frequency of allele A at marker M1 when Pa = 0.5 at 1% significance level. PD1 = 0.5; PD2 = 0.5; PB = 0.5; DAD1 = {min(PD1, PA) – PD1, PA}/2; DBD2 = 0.05; penetrance f = 0.5; the sample size of cases R = 1,000 and controls ...

Figure Figure44 illustrates the power against a linkage disequilibrium coefficient DAD1 with a fixed value of DBD2 = 0.05. The power is calculated under f = 0.5, PA = Pa = 0.5, PB = Pb = 0.5 and sample size 1,000 cases and 1,000 controls. The powers of ‘Dom [dbl union] Dom’ and ‘Threshold’ were identical and the power curves over the 6 models are very similar to those in figure figure22 and and3.3. Figure Figure55 shows the power against the disease allelic frequency pD1 at D*1 when pD2 is fixed pD2 = 0.5, which have the same parameters as figure figure44 except DAD1 = {min(PD1, PA) – PD1PA}/2. Likewise, the power of ‘Dom ∩ Dom’ is the highest and ‘Dom [dbl union] Rec’ outperforms ‘Threshold’. Surprisingly, the power of the ‘Modified’ model was the third highest power at pD1 < 0.32, but it leveled off at a power of approximately 0.6.

Fig. 4
Analytical power against DAD1 when DBD2 = 0.05 at 1% significance level. PD1 = 0.5; PD2 = 0.5; PA = 0.5; penetrance f = 0.5; the sample size of cases R = 1,000 and controls S = 1,000; the notations of models are the same as figure figure22.
Fig. 5
Analytical power against frequency of disease allele D1 when PD2 = 0.5 at 1% significance level. DAD1 = {min(PD1, PA) – PD1, PA}/2, DAD1 = DBD2 = 0.05; penetrance f = 0.5; the sample size of cases R = 1,000 and controls S = 1,000; the notations ...

Figure Figure66 shows the required sample size for cases and controls across the full range of penetrance f to achieve 80% power at 1% significance level when we assume the ratio of cases to controls is 1 (equal size). The sample size is calculated under PA = Pa = 0.5; PB = Pb = 0.5 and DAD1 = DBD2 = 0.05. As shown in figure figure6,6, at f > 0.5, all models need less than 1,500 case (1,500 control) subjects except the ‘Modified’, ‘Dom [dbl intersection] Dom’, and ‘Threshold’ model which need less than 1,000 case subjects. The ‘Modified’ model needs the largest sample size over the full range of f, followed by ‘Dom [dbl union] Rec’ and ‘Rec [dbl union] Rec’ models.

Fig. 6
The required sample size of cases against penetrance f for 80% power achievement at 1% significance level. PD1 = 0.5; PD2 = 0.5; PA = 0.5; PB = 0.5; DAD1 = DBD2 = 0.05; ratio of case to control is 1. Model notations are the same as figure figure ...

Power Comparisons

Logistic Regression versus Cochran-Armitage (CA) Regression

In order to compare the score test derived from the logistic regression with the F test and LRT from CA regression, we performed simulation studies using SNaP with two scenarios. The first scenario was defined by common variants of SNPs with MAFs at the two loci of 0.3 (PA = PB = 0.3), disease-associated allele frequencies of 0.2 (PD1 = PD2 = 0.2) for a more realistic situation. The LD between a disease locus and a marker is DAD1 = DBD2 = 0.6; R2AD1 = R2BD2 = 0.1. The second scenario was defined by rare variants of SNPs with a minor allelic frequency at the two loci of 0.05 (PA = PB = 0.05) and disease-associated allele frequencies of 0.5 (PD1 = PD2 = 0.5) that are increased due to the floor effect to have the similar LD, DAD1 = DBD2 = 0.6; R2AD1 = R2BD2 = 0.3. Both scenarios have a penetrance f = 0.2. For each model 2,000 datasets were simulated. Each dataset has 200 cases and 200 controls. Table Table55 shows the empirical power of the score test, F and LRT at the 1 and 5% significance level in the common variants; table table66 shows the results of the rare variants of SNPs. The power of the three tests in both scenarios is asymptotically equivalent across the 6 interaction models.

Table 5
Power (%) at 1% and 5% significance level at DAQ1 = DBQ2 = 0.6; R2AQ1 = R2BQ2 = 0.1 (PA = PB = 0.3) from 2,000 datasets, each has 200 cases and 200 controls
Table 6
Power (%) at 1% and 5% significance level at DAQ1 = DBQ2 = 0.6; R2AQ1 = R2BQ2 = 0.3 (PA = PB = 0.05) from 2,000 datasets, each has 200 cases and 200 controls

Allelic-Based Method versus Genotypic-Based Methods

We have compared the allelic-based approach with two genotypic-based approaches. One is a score test of a genotypic-based model using a genotypic score in a logistic regression framework as follows.

logit[p(yi=1|wiγ,Zi,k=(AA,BB),,Zi,k=(aa,bb))]=α+wiγ+k=18Zi,k=(jl,mn)βk,

where j, l = (A, or a), m, n = (B, or b) and

Zi,k=(jl,mn)={1ifGk=(jl,mn)0otherwise.

Under the null hypothesis of no interaction between two unlinked loci, H0 : βAABB = … = βaaBb = 0. The score test can be derived from the logistic regression, which is approximately equivalent with the F and LRT from the CA regression. The other genotypic-based method we have compared is multifactor dimensionality reduction (MDR). Due to the use of cross-validation to calculate a prediction accuracy, it is not directly comparable [17]. Therefore, we used a chi-square test of homogeneity for high risk and low risk groups between cases and controls and estimated empirical type I error rates from the 5,000 null datasets and the empirical power of each model from the 2,000 datasets after controlling for the type I error rates at each model in both scenarios [17]. Tables Tables55 and and66 presented the power of the score test for the allelic-based method and the genotypic-based methods. The tables illustrate that the allelic-based method outperformed both a genotypic-based method using a genotypic score and the MDR.

Application to Major Depression Disorder

We applied our method to Major Depression Disorder (MDD) case and control data provided by dbGaP [6]. The 1,754 cases with the disease endpoint of diagnosis and 1,763 controls without the MDD were genotyped on the Perlegen 600K platform. We selected 7 candidate genes related with a serotonergic system and a brain-derived neurotrophic factor (BDNF) associated with a pathway of MDD. There were 78 SNPs available from the 7 candidate genes (SLC6A4, HTR2A, TPH1, TPH2, ITGB3, COMT and BDNF). As quality control procedures, SNPs that have significant deviation from Hardy-Weinberg equilibrium (p value < 6.0 × 10−4) and a minor allele frequency less than 0.01 were removed. 64 tag SNPs using Haploview were selected to reduce the number of tests due to the issue of multiple comparison. The procedure for the proposed method to search for genes associated with the MDD through their interaction consists of multiple step. First, two genes among the 7 genes were selected, then we performed a two-way interaction analysis with an SNP from each gene. All SNPs in the two unlinked genes were tested in a pair-wise fashion. 1,470 SNP pair combinations from 21 two-gene combinations were tested. At each two-way interaction, we performed a score test, an F test and an LRT. Only significant SNP pairs that have p values less than 3.4 × 10−5 after adjustment for Bonferroni correction were compared with a model consisting only of a main effect of each SNP to distinguish real allelic interactions. The criteria to select the best model between an interaction model and a main effect model is (1) the p value of the global test is less than 3.4 × 10−5 and (2) the p value of the model comparison between the interaction model and the main effect model is less than 0.01, and (3) the smaller Akaike Information Criterion (AIC) from the following models.

Main effect:HA1:logit[p(yi=1|wi,Xi,A,Xi,B)]=μ+wiγ+Xi,AβA+Xi,BβB

Two-way:HA2:logit[p(yi=1|wi,Xi,AB,Xi,Ab,Xi,aB)]=μ+wiγ+Xi,ABβAB+Xi,AbβAb+Xi,aBβaB

Second, SNP pairs detected by an interaction at the first step were considered for three-way interaction analysis by adding an SNP from one of the 5 genes that were not in the two-way interaction model. This procedure is called a ‘forward selection’ procedure. Similarly to the two-way interaction analysis, we carried out a comparison test of the significant three-way interaction model with a main effects model. The best three-way interaction model was selected by the same three criteria. Finally, these steps were continued until no further high dimensional interaction models are detected.

Based on the 7 selected genes with 64 tagging SNPs, we identified two interacting SNPs: rs11030104 in BDNF gene and rs9526240 in HTR2A gene (p value = 2.3 × 10−5 for a score and LRT test; p value = 8.11 × 10−6 for a comparison test between HA2 of a two-way interaction model and HA1 of a main effects model). Two pairs of SNPs, rs6265 and rs11030104 in BDNF, and rs9526240 and rs6561335 in HTR2A are in LD with R2 = 0.89 and R2 = 0.8, respectively. These 4 pair-wise SNPs are also significantly interacting (p values are between 2.4 × 10−5 and 2.6 × 10−5). Based on the interacting SNPs, we performed a three-way interaction analysis and detected a three-way interaction of rs11030104 in BDNF, rs9526240 in HTR2A and rs652458 in TPH1 genes (p value for the interaction model is 5.0 × 10−5; p value for the comparison with a main effect is 0.002). The interaction of BDNF and HTR2A was identified in schizophrenia [18] and in unipolar depression [19]. The association of SLC6A4, HTR2A was detected in the MDD [20].

Discussion

In this study, a novel statistical approach to identify gene-gene interaction contributing to a disease trait in a case-control sampling design was proposed. In contrast with most of the currently available methods detecting interaction at the genotypic level, the proposed method searches for interaction at the allelic level with a new definition of interaction. Interaction occurs when the contribution to disease of a particular allele inherited in one gene depends on a particular allele inherited in other unlinked genes. Due to the new definition, the interpretation of result may be more straightforward at the allelic rather than at the genotypic level because the interacting alleles can be explained by combinations of disease-associated alleles. The score test derived from a logistic regression model is computationally efficient because one does not need maximum likelihood estimations. As shown in the simulation studies, the power of the score test is asymptotically equivalent to that of the F test and LRT derived from the CA regression. Furthermore, our approach allows an adjustment for non-genetic covariates as nuisance parameters, which may be necessary in complex trait analysis.

The terms ‘gene-gene interaction’ or ‘epistasis’ are used with different meanings within different areas of genetics [6]. Recently, Phillips reviewed the conceptual definition and the differences, and also pointed out that there are conceptual barriers to generating a more unified definition [6]. Epistasis was originally defined as a mutation (or an allele) that masks or stimulates the effects of the other mutation on a phenotype. Gene-gene interaction is a broad concept that encompasses a variety of genetic phenomena. In this paper, we use gene-gene interaction in a general way to involve the epistasis. Detection of the genetic interaction at the allelic level is conceptually similar to epistasis for which the biological interpretation is straightforward.

For the high-order interaction of multiple markers, the issue of multiple comparisons still remains either in independent SNPs or in dependent SNPs due to multiple procedures to distinguish the real allelic interaction from the genetic effects. As illustrated in the application with MDD, each procedure involves testing two or three hypotheses and many combinations of SNPs involve testing many hypothesis testings. Bonferroni correction for the proposed interaction models was used to find the significant sets of SNPs at each step depending on the number of SNPs being analyzed. For each comparison test, a threshold of 0.01 or less was used for the significance level. It has been known that the Bonferroni correction often produces conservative results when LD is ignored. Permutation analysis is an alternative approach to generate an empirical null distribution of each step; however, it may not be computationally feasible for real data applications because of the multiple hierarchical steps for high dimensional interactions. Some machine learning approaches use cross-validation to predict an accuracy rate and permutation analysis to estimate p value, both of which are computationally expensive.

Searching for gene-gene interaction is a necessary procedure in a genome-wide association study because a disease is a complicated network of multiple genes. Effort to develop computationally efficient methods to detect the interacting SNPs at the whole genome level is underway. Future work should extend this approach to the genomewide association study. Alternatively, an incorporation of biological knowledge of gene pathway to prioritize biologically important genes may reduce the search space and increase power to detect interacting genes associated with a disease.

Acknowledgments

We would like to thank two anonymous reviewers and Peter H. Baenziger for thoughtful critiques that have resulted in a much improved manuscript. This work was supported by NIH/NHLBI R01 HL095086-01. Funding support for Major Depression: Stage 1 Genomewide Association in Population-Based Samples was provided by NIH and the genotyping of samples was provided through the Genetic Association Information Network (GAIN). The dataset(s) used for the analyses described in this manuscript were obtained from the GAIN Database found at http://view.ncbi.nlm.nih.gov/dbgap-controlled through dbGaP accession number phs000020.v1.p1. Samples and associated phenotype data for Major Depression: Stage 1 Genomewide Association in Population-Based Samples were provided by P. Sullivan.

Appendix A

Let Uτ = (UAB, UAb, UaB)τ be a score function, which is the derivative of the log-likelihood function with respect to β = (βAB, βAb, βab) respectively

Uτ=Uβ(βH0,αˆ)=(UAB,UAb,UaB)τ=(logLβAB,logLβAb,logLβaB)τ=(xi,AB(yi-y¯),xi,Ab(yi-y¯)xi,aB(yi-y¯))τ

Let I (α, βAB, βAb, βaB) be the observed Fisher information matrix, which is the second derivative of the log-likelihood function with respect to (α, βAB, βAb, βaB) respectively

I(α,βAB,βAb,βaB)=y¯(1-y¯)(nxi,ABxi,Abxi,aBxi,ABxi,AB2xi,ABxi,Abxi,ABxi,aBxi,Abxi,ABxi,Abxi,Ab2xi,Abxi,aBxi,aBxi,ABxi,aBxi,Abxi,aBxi,aB2)=(IααIαβIβαIββ),

where Iαα = [partial differential]2 log L/[partial differential]α2, Iαβ = [partial differential]2 log L/[partial differential]α[partial differential], Iββ = [partial differential]2 log L/[partial differential]β2.

Appendix B

Under the assumption of no covariates for simplification, we calculate non-centrality parameter λscore = ΩτΣ−1Ω where Ω and Σ are the mean and covariance matrix of U* respectively. Let X = (XAB, XAb, XaB, Xab), where XAB = (X1,AB, X2,AB, …, Xn,AB)τ. yi = 1 if a subject is affected by a disease and yi = 0 if not affected. Let fklmn = P(y = 1 | D*1, D*2) be the joint penetrance for the two disease loci D*1 and D*2, k, l = (D1, d1), m, n = (D2, d2).

E(XABy)=klmnXAB×y×P(y=1,D1*=DkDl,D2*=DmDn,M1,M2)=klmnxAB×P(y=1|D1*,D2*,M1,M2)×P(D1*,D2*,M1,M2)=1×{f1111×P(D1*=D1D1,D2*=D2D2,M1=AA,M2=BB)++f2222×P(d1d1,d2d2,AA,BB)}+12×{f1111×P(D1D1,D2D2,AA,Bb)++f2222×P(d1d1,d2d2,AA,Bb)}+12×{f1111×P(D1D1,D2D2,Aa,BB)++f2222×P(d1d1,d2d2,Aa,BB)}+14×{f1111×P(D1D1,D2D2,Aa,Bb)++f2222×P(d1d1,d2d2,Aa,Bb)}=f1111{p(AD1)2P(BD2)2+12×2p(AD1)2P(BD2)P(bD2)+12×2p(AD1)p(aD1)P(BD2)2+14×4p(AD1)p(aD1)P(BD2)P(bD2)}++f2222{p(Ad1)2P(Bd2)2+12×2p(Ad1)2P(Bd2)P(bd2)+12×2p(Ad1)p(ad1)P(Bd2)2+14×4p(Ad1)p(ad1)P(Bd2)P(bd2)}

E(XABy)=P(AD1)P(BD2){f1111P(D1)P(D2)+f1112P(D1)P(d2)+f1211P(d1)P(D2)+f1212P(d1)P(d2)}+P(AD1)P(Bd2){f1112P(D1)P(D2)+f1122P(D1)P(d2)+f1212P(d1)P(d2)+f1222P(d1)P(d2)}+P(Ad1)P(BD2){f1211P(D1)P(D2)+f1212P(D1)P(d2)+f2211P(d1)P(D2)+f2212P(d1)P(d2)}+P(Ad1)P(Bd2){f1212P(D1)P(D2)+f1222P(D1)P(d2)+f2212P(d1)P(D2)+f2222P(d1)P(d2)}

Utilizing P (AD1) = DAD1 + PAPD1; P (Ad1) = DAd1 + PAPd1; P(BD2) = DBD2 + PBPD2; P(bd2) = Dbd2 + PbPD2 and fD1D1 = f1111P(D1)P(D2) + f1112P(D1)P(d2) + f1211P(d1) P (D2) + f1212P (d1) P (d2), we derived

E(XABy)=(DAD1+PAPD1)×(DBD2+PBPD2)fD1D2++(DAd1+PAPd1)×(DBd2+PBPd2)fd1d2=fPAPB+(fD1D2-fD1d2-fd1D2+fd1d2)DAD1DBD2+(fD1-fd1)DAD1PB+(fD2-fqd)DBD1PA.

E (XAby), E (XaBy) can be calculated in the same way.

Using E (XAB) = PAPB, E (XAb) = PAPb, E (XaB) = PaPB, E(Xab) = PaPb, E(X2AB) = PAPB (1 − 0.5 Pa)(1 − 0.5 Pb), E(XABXAb) = 0.5 PAPBPb (1 − 0.5 Pa), …, E (XaBXab) = 1/4 PaPBPb (1 − 2 PA),

E(xi,AB2xi,ABxi,Abxi,ABxi,aBxi,ABxi,Abxi,Ab2xi,Abxi,aBxi,ABxi,aBxi,Abxi,aBxi,aB2)=NRN(1-RN)(PAPB(1-0.5Pa)(1-0.5Pb)0.5PAPBPb(1-0.5Pa)0.5PAPBPb(1-0.5Pb)0.5PAPBPb(1-0.5Pa)PAPb(1-0.5Pa)(1-0.5PB)0.25PAPBPaPb0.5PAPBPa(1-0.5Pb)0.25PAPBPaPbPaPB(1-0.5PA)(1-0.5Pb))

Therefore, and can be expressed as follows.

Ω=E(U*)=E((iyixi,ABiyixi,Abiyixi,aB)-y¯(ixi,ABixi,Abixi,aB))τ=[RN(PAPBPAPbPaPB)+(fD1D2-fD1d2-fd1D2+fd1d2)(DAD1DBD2DAD1DbD2DaD1DBD2)+(fD1-fd1)(DAD1PBDAD1PbDaD1PB)+(fD2-fd2)(DBD2PADbD2PADBD2Pa)-RN(PAPBPAPbPaPB)]=[(fD1D2-fD1d2-fd1D2-fd1d2)(DAD1DBD2DAD1DbD2DaD1DBD2)+(fD1-fd1)(DAD1PBDAD1PbDaD1PB)+(fD2-fd2)(DBD2PADbD2PADBD2Pa)]

Σ=V(U*)=Iββ-IβαIαα-1Iαβ=y¯(1-y¯)(xi,AB2xi,ABxi,Abxi,ABxi,aBxi,ABxi,Abxi,Ab2xi,Abxi,aBxi,ABxi,aBxi,Abxi,aBxi,aB2)-1Ny¯2(1-y¯)2(xi,ABxi,Abxi,aB)(xi,ABxi,Abxi,aB)=NRN(1-RN)[(PAPB(1-0.5Pa)(1-0.5Pb)PAPBPb(1-0.5Pa)PAPBPa(1-0.5Pb)PAPBPb(1-0.5Pa)PAPb(1-0.5Pa)(1-0.5PB)PAPBPaPbPAPBPa(1-0.5Pb)PAPBPaPbPaPB(1-0.5PA)(1-0.5Pb))-RR(1-RN)(PAPBPAPbPaPB)(PAPBPAPbPaPB)].

Appendix C

We multiply both sides by X = (XAB, XAb, XaB, Xab), where XAB = (X1,AB, X2,AB, …, Xn,AB)τ and take the expectation in order to derive regression coefficients in equation (3)

E(XAByXAbyXaByXaby)=E(XAB2XABXAbXABXaBXABXabXAbXABXAb2XAbXaBXAbXabXaBXABXaBXAbXaB2XaBXabXabXABXabXAbXabXaBXab2)×(βABβAbβab)

The expectation of each element in matrix E (X′ X) is calculated with the assumption of no covariates as shown in Appendix A.

E(XABy)=klmnXAB×y×P(y=1,D1*=DkDl,D2*=DmDn,M1,M2)=klmnxAB×P(y=1|D1*,D2*,M1,M2)×P(D1*,D2*,M1,M2)=(DAD1+PAPD1)×(DBD2+PBPD2)fD1D2++(DAd1+PAPd1)×(DBd2+PBPd2)fd1d2=fPAPB+(fD1D2-fD1d2-fd1d2)DAD1DBD2+(fD1-fd1)DAD1PB+(fD2-fqd)DBD1PAβ=(αβABβAbβaB)=K-1[RN(1PAPBPAPbPaPB)+(fD1D2-fD1d2-fd1D2+fd1d2)(0DAD1DBD2DAD1DbD2DaD1DBD2)+(fD1-fd1)(0DAD1PBDAD1PbDaD1PB)+(fD2-fd2)(0DBD2PADbD2PADBD2Pa)],
(2)

where

K(1PAPBPAPbPaPBPAPBPAPB(1-0.5Pa)(1-0.5Pb)0.5PAPBPb(1-0.5Pa)0.5PAPBPa(1-0.5Pb)PAPb0.5PAPBPb(1-0.5Pa)PAPb(1-0.5Pa)(1-0.5PB)0.25PAPBPaPbPaPB0.5PAPBPa(1-0.5Pb)0.25PAPBPaPbPaPB(1-0.5PA)(1-0.5Pb)).

References

1. Ritchie MD, Hahn LW, Roodi N, Bailey R, Dupont WD, Parl FF, Moore J. Multifactor-dimensionality reduction reveals high-order interactions among estrogen-metabolism genes in sporadic breast cancer. Am J Hum Genet. 2001;69:138–147. [PubMed]
2. Cook NR, Zee RY, Ridker PM. Tree and spline based association analysis of gene × gene interaction models for ischemic stroke. Stat Med. 2004;23:1439–1453. [PubMed]
3. Chen SH, Sun J, Dimitrov L, Turner AR, Adams TS, Meyers DA, Chang BL, Zheng SL, Grönberg H, Xu J, Hsu FC. A support vector machine approach for detecting gene-gene interaction. Genet Epidemiol. 2008;32:152–167. [PubMed]
4. Kooperberg C, Ruczinski I. Identifying interacting SNPs using Monte Carlo logic regression. Genet Epidemiol. 2005;28:157–170. [PubMed]
5. Jung J, Sun B, Kwon D, Koller D, Foroud T: Allelic based gene-gene interaction association with quantitative trait loci. Genet Epidemiol, 2009, in press. [PMC free article] [PubMed]
6. Pillips PC. Epistasis – the essential role of gene interactions in the structure and evolution of genetic systems. Nat Rev. 2008;9:855–867. [PMC free article] [PubMed]
7. Boomsma DI, Willemsen G, Sullivan PF, Heutink P, Meijer P, Sondervan D, Kluft C, Smit G, Nolen WA, Zitman FG, Smit JH, Hoogendijk WJ, van Dyck R, de Geus EJ, Penninx BW. Genome-wide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects. Eur J Hum Genet. 2008;16:335–342. [PubMed]
8. Anderson JA. Separate sample logistic discrimination. Biometrika. 1972;59:19–35.
9. Prentice RL, Pyke R. Logistic disease incidence models and case-control studies. Biometrika. 1979;66:403–411.
10. Scott AJ, Wild CJ. Fitting regression models to case-control data by maximum likelihood. Biometrika. 1997;84:57–71.
11. Self SG, Mauritsen RH. Power/Sample size calculations for generalized linear models. Biometrics. 1988;44:79–86.
12. Falconer DS, Mackay TFC. Introduction to quantitative genetics. ed 4. London: Longman; 1996.
13. Armitage P. Tests for linear trends in proportions and frequencies. Biometrics. 1955;11:375–386.
14. Cochran WG. Some methods for strengthening the common tests. Biometrics. 1954;10:417–451.
15. McAleer M. Exact tests of a model against nonnested alternatives. Biometrika. 1983;70:285–288.
16. Nothnagel M. Simulation of LD block-structured SNP haplotype data and its use for the analysis of case-control data by supervised learning methods. Am J Hum Genet. 2002;71(suppl):A2363.
17. Millstein J, Conti DV, Gilliland FD, Gauderman WJ. A testing framework for identifying susceptibility genes in the presence of epistasis, Am J Hum Gene. 2006;78:15–27. [PubMed]
18. Alfimove MV, Alfimova MV, Lezheôko TV, Golimbet VE, Korovaôtseva GI, Lavrushkina OM, Kolesina NIU, Frolova LP, Muratova AA, Abramova LI, Kaleda VG. Investigation of association of the brain-derived neurotrophic factor (BDNF) and a serotonin receptor 2A (5-HTR2A) genes with voluntary and involuntary attention in schizophrenia. Zh Nevrol Psikhiatr Im S S Korsakova. 2008;108:62–69. [PubMed]
19. Kotte A, McQuaid JR, Kelsoe JR: Psychotherapeutic mechanisms of change: the role of genes in depression treatment outcome, Abstract, Am Soc Hum Genet, 2007.
20. Lazary J, Lazary A, Gonda X, Benko A, Molnar E, Juhasz G, Bagdy G. New evidence for the association of the serotonin transporter gene (SLC6A4) haplotypes, threatening life events, and depressive phenotype. Biol Psychiatry. 2008;64:498–504. [PubMed]

Articles from Human Heredity are provided here courtesy of Karger Publishers