Eigen-Epistasis for detecting gene-gene interactions

Abstract

Background

A large amount of research has been devoted to the detection and investigation of epistatic interactions in genome-wide association studies (GWASs). Most of the literature focuses on low-order interactions between single-nucleotide polymorphisms (SNPs) with significant main effects.

Results

In this paper we propose an original approach for detecting epistasis at the gene level, without systematically filtering on significant genes. We first compute interaction variables for each gene pair by finding its Eigen-Epistasis component, defined as the linear combination of Gene SNPs having the highest correlation with the phenotype. The selection of significant effects is done using a penalized regression method based on Group Lasso controlling the False Discovery Rate.

Conclusion

The method is tested against two recent alternative proposals from the literature using synthetic data, and shows good performances in different settings. We demonstrate the power of our approach by detecting new gene-gene interactions on three genome-wide association studies.

Background

Genome Wide Association Studies (GWASs) look for genetic markers linked to a phenotype of interest. Typically, hundreds of thousands of single nucleotide polymorphisms (SNPs) are studied for a limited number of individuals using high-density genotyping arrays. Usually the association between each SNP and the phenotype is tested using single-marker methods. Multiple markers may also be considered, but these are typically selected using simple forward-selection methods. GWASs are a powerful tool for investigating the genetic architecture of complex diseases and have been successful in identifying hundreds of variants. However, they have been able to explain only a small proportion of the phenotypic variations expected from classical family studies [1]. A number of explanations for this missing heritability have been put forward. For example, it has been suggested that shared environments among relatives are not adequately taken into account. Another suggestion is that much larger numbers of variants with small effects remain to be identified. Rare variants, which are difficult to find using existing genotyping arrays [1], seem to be important causal factors, and so do structural variations. But complex diseases may also be caused, at least in part, by complex genetic structures with multiple interactions between markers (a phenomenon termed epistasis). Whereas in pedigree studies the genetic effect on phenotype is seen as part of the additive genetic variance, in GWASs it is seen as an unmeasured interaction between genes [2]. For example, Zuk et al. proposed a model that takes into account epistatic interaction in relation to Crohn’s disease [3]. They found that 80% of the missing heritability could be due to genetic interactions.

In recent years a number of methods for studying epistasis have been proposed and reported in various reviews [4–6]. They vary in terms of their data analysis (genome-wide or filtering) and their statistical methodology (Bayesian, frequentist, machine learning or data mining). Most of them focus on single-locus interactions, but considering interactions at the gene level can have several advantages. First, given that genes are the functional unit of the genome, results may be more biologically interpretable. Second, genetic effects are more easily detected when SNP effects are aggregated together. Third, gene-based analysis simplifies the multiple testing problem by reducing the number of variables. Several gene-gene methods have been proposed. These are based on a summarizing step which is used to obtain information at the gene level. In more recent methods, filters or penalized models are used to make the method applicable to a large number of genes, while older methods are only applicable to two or a very limited number of genes. For the summarizing step, most methods resort to a principal components (PC) approach, but each method has its specific characteristics. We describe some of these below.

Chatterjee et al. harnessed Tukey’s one-degree-of-freedom method to investigate interaction between two genes [7]. Their method is based on the assumption that the SNPs included in each gene region act as surrogates for an underlying biological phenotype. The genotypic information for the gene region is extracted as a single component by a weighted sum of all SNPs. The weights are determined according to the SNP’s correlation with the trait. The product of the two sums is then introduced as the gene-gene interaction term into a logistic model, where marginal effects are represented by the respective sums. Building on this idea, Wang et al. compared two different interaction tests [8]. On the one hand, they used Principal Component Analysis (PCA) to summarize SNP information within a gene, and on the other hand they used Partial Least Squares (PLS) to extract components that summarize, first, the information among SNPs in a gene and, second, the correlation between SNPs and the outcome of interest. They then proposed an interaction test based on either the first PC or the first PLS component for each gene, and were able to show that the PCA and PLS methods often outperformed Tukey’s one-degree-of-freedom method. But it is worth noting that the main objective of these three methods was improving the detection of associations in the presence of gene-gene interactions, rather than identifying the interactions themselves. Other approaches based on principal component analysis have since been proposed for epistasis detection. Li et al. proposed selecting, as the gene representation, PCs that are able to explain at least 80% of the variation [9].

Genotypic data are characterized by the high correlation among markers resulting from so-called linkage disequilibrium (LD). Procedures that take LD information into account have been developed for epistasis detection. For example, He et al. proposed an approach using LD information to weight genotype scores which are then aggregated using principal components [10]. Rajapakse et al. developed a gene-based test of interactions for case-control studies which compares LD patterns between cases and controls [11]. Using the same idea, Peng et al. used a canonical correlation-based U-statistic model (CCU) to detect co-association in case-control studies [12]. The idea is to test for two given genes the difference between canonical correlation coefficient computed by Canonical Correlation Analysis (CCA) among cases and among controls. Their work was subsequently extended to include kernel [13, 14].

However most of these methods can be applied only to a reduced number of genes. Computational constraints mean that it is not feasible to model all gene-gene interactions directly. One way of overcoming this is to reduce the gene-gene search space by eliminating unimportant genes, and to this end two-step procedures have been developed that first filter out specific genes or SNPs through a genome-wide search before testing for interactions. One example of this is the model-based kernel machine method (3G-SPA) proposed by Li and Cui, which first performs a search for gene pairs contributing to the overall phenotypic variations [15]. Significant pairs are then tested for interaction effects. Another attractive alternative is offered by penalized regression methods that select a subset of important predictors out of a large number of potential predictors. These methods operate by shrinking the size of the coefficients. The coefficients of predictors with little or no apparent effect are force to be set to zero, reducing the effective degrees of freedom and in many cases making model selection possible. A few approaches using penalized models have been proposed. D’Angelo et al. combined principal component analysis and lasso penalized regression [16]. Wang et al. used a principal component analysis combined with an L1 penalty, with adaptive weights based on gene size, pathway support and effect size [17].

Here we propose a Group Lasso approach [18] that takes into account the group structure of each gene in order to detect epistasis. We introduce Gene-Gene Eigen-Epistasis (G-GEE) as a new approach for computing the gene-gene interaction part of the model, and we compare G-GEE with two different interaction variable modeling approaches inspired by previous proposals in the literature, namely PCA and PLS. An adaptive ridge-cleaning approach [19] is then used in order to compute p-values for each group.

In the next section, we detail each model and outline the design of the simulation studies performed to compare the performance of the different approaches. In the Results section, the findings of the simulation studies are shown, and we illustrate our approach on three real datasets relating to ankylosing spondylitis, thyroid carcinomas and inflammatory bowel disease. The different approaches and the results are discussed in the last section.

Methods

We consider n individuals where y=(y ₁,y ₂,…,y _n )^T denotes the vector of trait values. For each individual, genetic variants among G genes are considered. Each gene is described by a given number of SNPs p _g where $\sum _{g} p_{g}=p$. The SNP matrix $\mathbf {X} \in \mathbb {R}^{n \times p}$ considers an additive coding scheme in which the genotype value of each SNP j from individual i is denoted X _ij ∈{1,2,3}. X _i is a p-dimensional vector of covariates for observation i and for j∈{1,…,p}. X ^g denotes the submatrix of X whose columns are the p _g SNPs of gene g. A generalized linear model is generally assumed for GWAS, where the phenotype is considered as a random variable y _i whose conditional expectation can be written as a function of the covariates X _i and their interactions Z _i ,

$$g\left(E \left[y_{i} | X \right]\right) = \boldsymbol{X}_{i}^{T} \boldsymbol{\beta} + \boldsymbol{Z}_{i}^{T} \boldsymbol{\gamma}, $$

where

$$\boldsymbol{\beta} = \left(\underbrace{\beta_{1,1}, \beta_{1,2},\ldots, \beta_{1,p_{1}}}_{gene_{1}},\ldots, \underbrace{\beta_{G,1},\ldots, \beta_{G,p_{G}}}_{gene_{G}} \right)^{T}, $$

and Z _i is the ith line of the matrix of interactions and γ a parameter vector of appropriate dimension. When the phenotype is binary (case control study), it is usual to assume a logistic model where g() is the logit and Y is assumed to follow a binomial distribution. Below we will consider only quantitative phenotypes using a classical linear model. In this case g() is the identity and the residuals are assumed to be Gaussian.

The main effect of each gene is modeled through the sum of the effects of all its SNPs. Concerning interaction effects, we compute new variables representing interactions between two specific genes and define as a group all the interaction variables related to a given pair of genes. The matrix of interaction is thus structured into G(G−1)/2 submatrices:

$$\boldsymbol{Z}= \left[\boldsymbol{Z}^{11} \cdots \boldsymbol{Z}^{rs} \cdots \boldsymbol{Z}^{G(G-1)/2}\right] $$

where Z ^rs describes the interactions between the two genes r and s. The parameter vector γ is accordingly structured into sub-vectors γ ^rs. We will now present and compare three different approaches for modeling gene-gene interactions.

Modeling gene-gene interactions

Let us consider two genes r and s described respectively by p _r and p _s SNPs. A possible interaction term describing the epistasis between the two genes is

$$ {\boldsymbol{Z}_{i}^{rs}}^{T} \boldsymbol{\gamma}^{rs} = \sum_{j=1}^{p_{r}} \sum_{k=1}^{p_{s}} \gamma_{jk}^{rs} X_{ij}^{r} X_{ik}^{s}. $$

(1)

We hereafter set $\boldsymbol {W}^{rs}=\{ X_{ij}^{r} X_{ik}^{s} \}_{i=1\cdots n}^{j=1,\cdots,p_{r}; k=1,\cdots,p_{s}}$. In this case the submatrix of interactions is Z ^rs=W ^rs and $\boldsymbol {\gamma }^{rs}=\{\gamma _{jk}^{rs}\}$ is a vector of size p _r p _s . The number of parameters in such a model is obviously too large to be reliably estimated. For this reason a number of papers in the literature consider reducing the dimension of γ.

In this paper we will consider three different methods for reducing the dimension reduction, namely Principal Component Analysis (PCA), Partial Least Squares (PLS), and our proposed Gene-Gene Eigen-Epistasis approach that we have termed G-GEE.

Principal component analysis

Principal Component Analysis (PCA) can reduce the number of variables describing each gene r from p _r to q _r <p _r . Considering gene r described by p _r SNPs, we compute the matrix of the first q principal components

$$C^{r} = \boldsymbol{X}^{r} U^{r}, $$

where U ^r is the matrix of the first q _r principal axes. Using C ^r and C ^s instead of X ^r and X ^s in the computation of the interaction allows the number of parameters relative to each interaction to be controlled. This control is achieved by choosing the number of principal components q. The PCA model that we describe draws upon ideas in [20]. The interaction term takes the form

$${\boldsymbol{Z}_{i}^{rs}}^{T} \boldsymbol{\gamma}^{rs} = \sum_{j=1}^{q} \sum_{k=1}^{q} \gamma_{jk}^{rs} C_{ij}^{r} C_{ik}^{s}. $$

Relating this expression to the general form of the interaction term $\boldsymbol {W}_{i}^{rs}$ described above, we can see that performing PCA prior to computing the interactions is a means of constraining the linear interaction term of Eq. 1.

The submatrix of interactions is $\boldsymbol {Z}^{rs}=\{ C_{ij}^{r} C_{ik}^{s} \}_{i=1\cdots n}^{j=1, \cdots, q;k=1, \cdots, q }$, and $\boldsymbol {\gamma }^{rs}=\{\gamma ^{rs}_{jk}\}$ is a vector of size q ² describing the interaction between genes r and s. In particular, if a single principal component is chosen, there will be only one parameter to estimate per interaction.

Partial least squares

Wang et al. proposed an alternative method for integrating interactions using a PLS approach [8]. Let (X ^r,X ^s) be the genotypic matrix for the given pair of genes (r,s). Their approach computes the components that maximize cov ²(X ^r u,T v), with T=(y,X ^s) and (u,v) the weight vectors. The interaction of a couple of genes (r, s) is then represented by the first q components:

$$ {\boldsymbol{Z}_{i}^{rs}}^{T} \boldsymbol{\gamma}^{rs} =\sum^{q}_{j=1} \gamma^{rs}_{j} T^{rs}_{ij}. $$

In this approach phenotypic information is retained when the interaction variables are constructed.

Gene-gene Eigen-Epistasis

We propose an original approach for modeling interactions. The general idea is to consider the interaction variable between the two genes r and s as a function f _u (X ^r,X ^s) parameterized by u. One way to estimate u is to maximize the correlation between the interaction function and the phenotype:

$$\hat{\boldsymbol{u}}= arg \max_{\boldsymbol{u}, \| \boldsymbol{u}\|=1} cov^{2}\left(\boldsymbol{y},f_{\boldsymbol{u}}\left(\boldsymbol{X}^{r},\boldsymbol{X}^{s}\right)\right). $$

If we consider the function f to be linear, our problem becomes easily tractable and has only one solution. Setting

$$ \boldsymbol{Z}^{rs}=f_{\boldsymbol{u}}\left(\boldsymbol{X}^{r},\boldsymbol{X}^{s}\right) = \boldsymbol{W}^{rs}\boldsymbol{u}, $$

where $\boldsymbol {W}^{rs}=\{ X_{ij}^{r} X_{ik}^{s} \}_{i=1\cdots n}^{j=1 \cdots, p_{r}; k=1,\cdots, p_{s}}$ and $\boldsymbol {u} \in \mathbb {R}^{p_{r} p_{s}}$ we obtain the following problem:

$$ \begin{aligned} \max_{\boldsymbol{u},\|\boldsymbol{u}\|=1}\left\|\hat{cov}\left[{\boldsymbol{W}^{rs}}\boldsymbol{u},\boldsymbol{y}\right]\right\|^{2}&= \max_{\boldsymbol{u},\| \boldsymbol{u}\|=1}\left\|\boldsymbol{u}^{T} {\boldsymbol{W}^{rs}}^{T} \boldsymbol{y}\right\|^{2}\\ &=\max_{\boldsymbol{u},\| \boldsymbol{u}\|=1}\boldsymbol{u}^{T} {\boldsymbol{W}^{rs}}^{T}\boldsymbol{y}\boldsymbol{y}^{T} \boldsymbol{W}^{rs} \boldsymbol{u}. \end{aligned} $$

(2)

The solution u is the eigenvector corresponding to the largest eigenvalue of the matrix W ^{rs T} y y ^T W ^rs, which is the vector W ^{rs T} y. The complexity of computing u is therefore in O(np _r p _s ). We then use the projection of the matrix W ^rs on u as the interaction variable. The resulting Eigen-Epistasis vector Z is the linear combination of all the SNP-SNP interactions being the most correlated with the phenotype. In its construction, G-GEE has similarities with PLS. The main difference lies in the original design matrix. PLS searches for components that maximize cov ²(X ^r u,y X ^s v), whereas G-GEE retains the component that maximizes cov ²(y,W ^rs u), with W ^rs the matrix of all pairwise interaction between the two genes r and s. Like PLS, G-GEE takes phenotypic information into account in the construction of the interaction variables. Other methods as such as CCU [12] and the kernel versions of CCU [13, 14] that we referred to in the introduction also consider the phenotype in their construction, but these methods can be applied only to case-control problems.

Estimation of coefficients

We propose a Group Lasso approach [18] for estimating the parameters of linear or logistic (case control) regression. A group comprises either the SNPs of a given gene, or interaction terms relative to a given gene-pair interaction. In the particular case of linear regression, the model parameters are estimated by:

$$\begin{aligned} \hat{\boldsymbol{\theta}}&=\left(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}}\right) = {\underset{\boldsymbol{\beta},\boldsymbol{\gamma}}{\arg\!\min}}\left(\sum_{i}{\left(y_{i} - \boldsymbol{X}_{i} \boldsymbol{\beta} - \boldsymbol{Z}_{i} \boldsymbol{\gamma}\right)^{2}} +\lambda \left[ \sum_{g} \sqrt{p_{g}}\|\boldsymbol{\beta}^{g}\|_{2}\right.\right.\\ &\qquad\qquad\quad\qquad\left.\left.+ \sum_{rs} \sqrt{p_{r} p_{s}}\|\boldsymbol{\gamma}^{rs}\|_{2}\right] \right), \end{aligned} $$

The parameter λ is selected by cross-validation.

In order to improve estimation accuracy and to obtain p-values for each of the selected groups, we use the adaptive ridge cleaning approach proposed by Bécu et al. [19]. This screen and clean procedure is a two-stage method. The group lasso model is first fitted on half of the data. The coefficient of the candidate groups selected by the model are then introduced into a ridge regression model fitted on the second half of the data with a specific penalty that allows the group structure to be taken into account. For each group the significance of the regression coefficients is estimated using permutation tests.

Simulation design

To evaluate the performance of the proposed approach, we conducted two simulation studies, the first using simulated data and the second using a real dataset relating to ankylosing spondylitis. In each case we compared the proposed G-GEE model to the two other interaction variable modeling approaches. The first simulation corresponds to a simplified context where all parameters were controlled and external interference limited, while the second simulation corresponds to a realistic context with a realistic pattern of minor allele frequency (MAF) and LD.

Design

Genotypes Our first (simplified) simulation study was adapted from the model used in [21] with an extension to control the MAF of each SNP. The n lines of the genotype matrix are an i.i.d. sample from a multivariate random vector $\boldsymbol {X}_{i}\sim \mathcal {N}_{p} (\boldsymbol {0},\boldsymbol {\Sigma })$. The correlation matrix Σ is block diagonal, each block corresponding to a gene. Two variables belonging to the same gene are correlated at level ρ=0.8 while all other correlations are null. Each SNP (column of the genotype matrix) is randomly assigned an MAF p from a uniform distribution between 0.05 and 0.5. An MAF value of 0.2 is assigned to all causal SNPs. The genotype frequencies derived from the Hardy-Weinberg equation are then used to discretize X _ik values to 0, 1 or 2. In practice, X _ik is set to 1 if $X_{ik}<q_{p^{2};N(0,1)}\phantom {\dot {i}\!}$, X _ik is set to 3 if $X_{ik}<q_{(1-p)^{2};N(0,1)}\phantom {\dot {i}\!}$ and X _ik is set to 2 otherwise.

In the second (realistic) simulation study using a real ankylosing spondylitis dataset, genes are randomly selected. The number of SNPs composing each genes varies according to the selection.

Phenotypes For both simulation studies, we generated phenotype vectors using two different schemes. Our first scheme corresponds to the model proposed by Wang et al. [17] (which, for the sake of brevity, we will refer to hereafter as the “Wang Pathway” model):

$$ \begin{aligned} Y_{i} = \beta_{0} + \sum_{g} \beta_{g} \left(\sum_{k \in \mathcal{C}} X_{ik}^{g} \right)+ \sum_{rs} \gamma_{rs} \left(\sum_{(j,k) \in \mathcal{C}^{2}} X_{ij}^{r} X_{ik}^{s} \right) + \epsilon_{i}, \end{aligned} $$

(3)

where $\mathcal {C}$ and $\mathcal {C}^{2}$ are respectively the set of causal SNPs and causal interactions, and ε _i a random Gaussian variable. For each causal gene g, we consider two causal SNPs and a coefficient β _g is assigned to the standardized sum of these causal SNPs. In the same way, for the interactions, all the causal SNPs from a causal pair (r,s) are pairwise multiplied and a coefficient γ _rs is assigned to the standardized sum of the product.

Our second scheme for simulating phenotypes is based on the following model:

$$ Y_{i} = \beta_{0} + \sum_{g} \beta_{g} \left(\sum_{k \in \mathcal{C}} X_{ik}^{g} \right) + \sum_{rs} \gamma_{rs} \left(\sum_{(j,k) \in \mathcal{C}^{2}} C_{ij}^{r} C_{ik}^{s} \right) + \epsilon_{i}. $$

(4)

The difference with the first model concerns the simulation of the interaction effect. In the second model the interaction effect for a causal couple (r,s) is defined as the product of the first PCA component $\boldsymbol {C}_{.1}^{r}$ of gene r and the first PCA component $\boldsymbol {C}_{.1}^{s}$ of gene s.

In both models, β ₀ is set to 0, and ε _i are generated independently from a $\mathcal {N}(0,\sigma ^{2})$, with σ ² determined from the coefficient of determination R ² that calibrates the strength of the association. Both simulation models can be written as $y_{i} = \boldsymbol {X}_{i}^{T} \boldsymbol {\beta } + \boldsymbol {Z}_{i}^{T} \boldsymbol {\gamma } + \epsilon _{i}$ where X the marginal effect genotype matrix and Z the interaction effect matrix.

Let us denote $\boldsymbol {Q} \boldsymbol {\phi }= [\boldsymbol {X}, \boldsymbol {Z}] \left [\begin {array}{c} \boldsymbol {\beta } \\ \boldsymbol {\gamma } \end {array}\right ]$ and

$$\begin{aligned} R^{2}&=\frac{\sum\left(\boldsymbol{Q_{i}\phi} - \bar{y}\right)^{2}}{\sum\left(\boldsymbol{Q_{i}\phi} + \epsilon_{i} - \bar{y}\right)^{2}}\\ &=\frac{\sum\left(\boldsymbol{Q_{i}\phi} - \bar{y}\right)^{2}}{\sum\left(\boldsymbol{Q_{i}\phi} - \bar{y}\right)^{2} + \sum\epsilon_{i}^{2} + \sum2\left(\epsilon_{i}\left(\boldsymbol{Q_{i}\phi}-\bar{y}\right)\right)}\\ &=\frac{\sum\left(\boldsymbol{Q_{i}\phi} - \bar{y}\right)^{2}}{\sum\left(\boldsymbol{Q_{i}\phi} - \bar{y}\right)^{2} + \textrm{n }\hat{\text{var}}(\epsilon_{i}) + 2 \mathrm{n} \hat{\text{cov}}\left(\epsilon_{i},\boldsymbol{Q_{i}\phi}-\bar{y}\right)}. \end{aligned} $$

We remark that:

$$\begin{aligned} 2 \mathrm{n} {\text{cov}}\left(\epsilon_{i},\boldsymbol{Q_{i}\phi}-\bar{y}\right) &= 2\mathrm{n cov}\left(\epsilon_{i},\boldsymbol{Q_{i}\phi}-\frac{\sum_{j} y_{j}}{n}\right)\\ &=2\mathrm{n cov}\left(\epsilon_{i},\boldsymbol{Q_{i}\phi}\right)-\sum_{j} \frac{2\mathrm{n}}{\mathrm{n}}\text{cov}\left(\epsilon_{i},y_{j}\right)\\ &=0 - 2\text{cov}\left(\epsilon_{i},\epsilon_{i}\right)=-2\sigma^{2} \end{aligned} $$

Thus, replacing $\hat {\text {var}}(\epsilon _{i})$ by σ ², and $\hat {\text {cov}}(\epsilon _{i},\boldsymbol {Q_{i}\phi }-\bar {y})$ by −σ ²/n, we obtain $R^{2} \approx \frac {\sum (\boldsymbol {Q_{i}\phi } - \bar {y})^{2}}{\sum (\boldsymbol {Q_{i}\phi } - \bar {y})^{2} + \mathrm {n}\sigma ^{2} -2\sigma ^{2}}$. This relation between R ² and σ ² gives us an expression for σ ² that depends on R ², $\sigma ^{2} = \frac {(R^{2} - 1) \sum (\boldsymbol {Q_{i}\phi } - \bar {y})^{2}}{R^{2} (2-n)}$.

We looked at how much of the coefficient of determination R ² is explained by main effects, and how much is explained by interaction effects, in order to determine their respective roles in the model.

For a similar reason, when simulating phenotypes, Wang et al. [17] examined how much of partial R ² was due to interaction effects. They selected coefficient values so that 30% of the partial R ² was explained by interaction effects. Li and Cui [15] did not use the R ² directly, but they simulated data assuming different proportions of interaction effects among the total genetic variance. In our study, once the phenotype y had been set for each simulated design matrix, we computed how much of the R ² could be attributed to interaction and main effects as $p_{I}=\frac {R^{2}_{I}}{R^{2}_{T}}$ and $p_{M}=\frac {R^{2}_{M}}{R^{2}_{T}}$ respectively, with $R^{2}_{I}$ the R-square value for the model containing only simulated interaction effects, $R^{2}_{M}$ the R-square value where there were only simulated main effects, and $R^{2}_{T}$ R-square value where there were both simulated main effects and simulated interaction effects.

Scenarios In the first (simplified) simulation study, genotypes are simulated as described in the design. We considered six genes, each composed of six SNPs and for 600 subjects. We define one causal interaction between genes and two causal genes with main effects, and the simulation takes place using two alternative simulation settings:

• (1) one interaction and two main effects involving the same genes

• (2) one interaction and two different main effects

For these two settings, different coefficients of determination, from 0.05 to 0.7, are considered and 1000 iterations are performed.

In the second (realistic) simulation study, genotypes come from a real dataset comprising 763 individuals. At each iteration we randomly select six genes of various size (from 1 to 1119 with a median of 2 SNPs) in the dataset. We consider the five following settings:

• (1) one interaction and two main effects involving the same genes

• (2) one interaction and two different main effects

• (3) one interaction effect only

• (4) two main effects only

• (5) no effects

For each setting, coefficients of determination, from 0.1 to 0.4, are considered and 500 iterations are performed.

For both simulation studies, main effects and interaction effects are weighted with the same coefficient values (β _g =γ _rs =2,∀g,r,s). For each interaction, the power is estimated as the proportion of detected interactions over the total number of simulations.

Real data illustration

To illustrate our approach we applied the proposed method on three real datasets related to ankylosing spondylitis, thyroid carcinomas and inflammatory bowel disease.

The dataset regarding ankylosing spondylitis consists of the French subset of the large study of the International Genetics of Ankylosing Spondylitis (IGAS) study [22]. For this subset, unrelated cases were recruited through the Rheumatology clinic of Ambroise Paré Hospital (Boulogne-Billancourt, France) or through the national self-help patients’ association: "Association Française des Spondylarthritiques". Population-matched unrelated controls were obtained from the "Centre d’Etude du Polymorphisme Humain", or were recruited as healthy spouses of cases. The protocol was reviewed and approved by the Ethics committee of the Ambroise Paré hospital. All participants gave their informed consent to the study. The application on thyroid carcinomas was carried out on a public dataset that came from the study of Luzón-Toro et al. on identification of epistatic interactions in two different types of thyroid carcinomas [23]. Finally, we used the Wellcome Trust Case-Control Consortium genome-wide association dataset to study Inflammatory Bowel Disease.

Results

Simulation studies

In the following, we will refer to the different simulation settings by using letters as described in Table 1.

Table 1 Effects simulated in each settings and referring names according to the phenotype simulation model

Results from the simplified simulation study

Figure 1 shows results obtained for the two settings. The first column gives the estimated power to detect the gene interaction as a function of the R ² values. The last two columns show heatmap matrices reflecting the proportion of significant values for each variable and each method over the 1000 simulations for different R ² values.

Fig. 1

Power and discoveries under a simplified context. The figures in the first column shows the power to detect interaction effects of the three methods depending on the R ². The last two columns show the ratio of the number of times where each variable was significant to the total number of simulations for a given R ². The panels a, b, c and d refer to the different simulation settings described in Table 1

In the first setting (Fig. 1 a, b), we consider genes 1 and 2, both having main and interaction effects. When the phenotype is simulated using the Wang Pathway model, the G-GEE and PLS methods have a higher power to detect the interaction effect than PCA method, which tends to identify only the two main effects of the two genes (Fig. 1 a). Whereas for PCA and PLS the power is nondecreasing with R ², for G-GEE we observe a U-shaped curve. For the smallest R ² values, which correspond to the most difficult cases, the power of G-GEE to detect the interaction tends to decrease. When R ² values reach 0.4, G-GEE’s power to detect the interaction starts to increase. The situation is different for the main effects, since G-GEE’s power to detect these increases continuously with R ² [see Additional file 1]. For PLS, the power to detect the interaction effect is continuously nondecreasing. Note, however, that for this method one of the two main effects (here gene 1) is detected to the detriment of the second, regardless of the value of R ². In the PCA phenotype simulation model (Fig. 1 b), G-GEE has a higher power than the other methods to detect interaction effects while retaining a good specificity, whatever the value of R ². The reasonably high power of the PCA method can be explained by the similarity between the phenotype simulation model and the estimation model. It is worth noting that in this first setting, only a few variables are falsely significant, which reflects a good specificity for all methods (the worst being for the gene 3 × gene 4 interaction variable in the case of the Wang Pathway model and r ²=0.1, where the false discovery rate is 0.068).

In the second setting (Fig. 1 c, d), genes 1 and 2 have only main effects, and genes 3 and 4 have only an interaction effect. When the phenotype is simulated using the Wang Pathway model, the interaction power of G-GEE is uniformly higher than that of the other methods (Fig. 1 c). For all values of R ², PCA tends to detect false main effects for genes 3 and 4, but not to detect interaction effects. In the PCA phenotype simulation (Fig. 1 d), PCA has a good power to detect interaction effects, but once again these good performances can be explained by the similarity between the simulation model and the estimation model. The interaction power for G-GEE is lower, but still good. With this model, only G-GEE tends to attribute a false main effect to genes 3 and 4. In this second setting, whatever the phenotype simulation model, the power of the PLS method is almost null. PLS identifies only the first gene as having a main effect, while the effects of genes 3 and 4 are not detected, whether as main or as interaction effects. Moreover, PLS tends to attribute a false interaction effect between genes 1 and 2.

To evaluate the performances of the different methods in a more complex context, we also consider a setting where we simulate 25 genes with four causal interactions between genes, and two genes with causal main effects. In these simulations, interaction genes are different from main effect genes, and we only consider the case where R ²=0.7. The results of this setting reflect the good performance of the G-GEE method over PCA and PLS in detecting interaction in a context where further interactions and different main effects are simulated [see Additional file 2].

Results from the realistic simulation study

Figures 2 and 3 show results for the first three settings. Fig. 2 shows the power to detect gene interaction depending on the R ², and Fig. 3 shows heatmaps of significant effects when R ²=0.2. In both figures, the upper row relates to phenotypes simulated using the Wang Pathway model, and the lower row to phenotypes simulated using the PCA model.

Fig. 2

Power under a realistic context. The figures show the power to detect interaction effects of the three methods depending on R ². The panels a, b, c, d, e and f refer to the different simulation settings described in Table 1

Fig. 3

Discoveries under a realistic context. Heatmaps of the ratio of the number of times where each variable was significant to the total number of simulations for R ²=0.2. The panels a, b, c, d, e and f refer to the different simulation settings described in Table 1

In the first setting (Fig. 2 a, b, Fig. 3a, b), the same two genes (genes 1 and 2) are simulated with main and interaction effects. In this setting G-GEE has the best power to detect interaction effects for all R ² values. The interaction power of PLS remains close to 0.3. The power of PCA depends on the phenotype simulation model. When the phenotype is simulated using the Wang Pathway model, the power is similar to PLS. When it is simulated using PCA it increases continuously, because of the similarity between the phenotype simulation model and the estimation model. Looking at the heatmaps (Fig. 3a, b) we can see that only a few variables are falsely significant. We also observe that unlike G-GEE, PCA and PLS can detect the two simulated main effects, with a preference for gene 1 in the case of PLS. These results are obtained when R ²=0.2, but other R ² values give similar results [see Additional file 3].

In the second setting (Fig. 2 c, d, Fig. 3c, d), main effects are simulated for genes 1 and 2, and one interaction is simulated between genes 3 and 4. In this setting the power of PLS to detect interaction effects is almost null, while the respective powers of PCA and G-GEE are different, according to which phenotype simulation model is used. Both methods have a higher power when the phenotype is simulated using the PCA model. Regarding the detection of main effects, the results are similar to the first setting, with G-GEE less successful than PCA and PLS (Fig. 3c, d). But unlike in the first setting, here some variables are falsely significant. False detections among interaction variables are more pronounced for G-GEE and concern genes that have been simulated to have only main effects. False detections among main effects are more pronounced for PCA and PLS when the phenotype is simulated using the Wang Pathway model and concern genes that have been simulated to have an interaction effect. Under the PCA phenotype model, false detections among main effects are more pronounced for PLS and G-GEE when R ² values are higher [see Additional file 3].

In the third setting, where only one interaction is simulated between genes 1 and 2, G-GEE has a higher power to detect interaction than PLS and PCA when the phenotype is simulated using the Wang Pathway model (Fig. 2 e). The power of PCA is higher in the PCA phenotype simulation model because of its similarity to the estimation model, whereas the power of PLS is almost null (Fig. 2 f). In the Wang Pathway phenotype simulation model, PCA and PLS both falsely detect main effects. In the PCA phenotype simulation model, the false detections are made by PLS and G-GEE (Fig. 3e, f). In all cases these false detections concern genes that are simulated to have an interaction effect.

Figure 4 shows the results for the fourth and fifth settings. The heatmap on the left corresponds to the fourth setting, where only two main effects are simulated. We remark that all methods successfully identify the main effect, PCA and PLS doing so with a higher power. False detections corresponding to the respective interaction effects are observed for G-GEE, and to a lesser extent for PLS. The figure on the right corresponds to the fifth setting, where no specific effects are simulated and the result shows that all three methods perform well with very few false detections.

Fig. 4

Discoveries for the fifth and the sixth settings. Heatmaps of the ratio of the number of times where each variable was significant to the total number of simulations for R ²=0.4 when only main effects are simulated for gene 1 and gene 2 (left), and when no effects are simulated (right)

In all settings, estimating the coefficients with the group lasso is more computationally expensive than constructing the interaction variables. G-GEE and PCA are quite similar in terms of computation time, whereas in some settings PLS has a slightly greater execution time than other methods. Note that the time required by G-GEE for constructing the interaction variables varies according to the number of SNPs that constitute each gene [see Additional file 4].

Percentage of R ² attributable to interaction and main effects respectively

Using each setting in both simulation studies, we determine the p _I and p _M average values that correspond to the proportion of the R ² attributed to interaction and main effects, respectively. For most settings, the p _I depends on the number of simulated effects. With one interaction and two main effects simulated the R ² part attributable to interaction effects is around 33% (Table 2 (B, C, D), Table 3 (C, D)). For the setting with numerous effects [see Additional file 2], the average p _I is 67% because we consider four interaction effects for only two main effects. Finally, as expected, when only interaction effects are simulated, the average p _I is 100% (Table 3 (E, F)) and 0% when only main effects are simulated (Table 3 (OME)). However, the R ² distribution between main and interaction effects is not distinguishable in the setting where the phenotype is simulated using the Wang Pathway model with the same main and interaction effects. The p _I and p _M values are all above 90% (Table 2 (A), Table 3 (A)). In the second simulation study, the R ² distribution is also not well divided between main and interaction effects when the phenotype is simulated under the PCA model, though p _M is still higher than p _I (Table 3 (B)).

Table 2 Average percentage of R ² attributable to interaction and main effects, by setting, in the first simulation study

Table 3 Average percentage of R ² attributable to interaction and main effects, by setting, in the second simulation study

Real data illustrations

Ankylosing spondylitis

Ankylosing spondylitis (AS) is a common form of inflammatory arthritis predominantly affecting the spine and pelvis. It occurs with a prevalence of 0.1% to 1.4% depending on the considered population [24]. Genetic factors account for more than 90% of the risk of susceptibility to AS. Human leukocyte antigen (HLA) class I molecule HLA B27, belonging to the Major Histocompatibility Complex (MHC) region, was the first genetic risk factor identified as associated with ankylosing spondylitis in the 1970’s [25, 26] and remains the most important risk locus for this pathology. Despite the strong association only a small portion of HLA-B27 carriers develop the disease. Furthermore, studies in families suggest that less than 50% of the overall genetic risk is due to HLA-B27, which suggests that other genetic factors are involved [27]. A number of updated reviews on AS genetics, including genome-wide association study (GWAS) results, identified new ankylosing spondylitis-susceptibility genes outside of the MHC region [28, 29].

We applied all the methods described above to the AS dataset. The data contain 408 cases and 358 controls, and each individual was genotyped for 116, 513 SNPs with Immunochip technology. For each SNP we obtained detailed genetic information, such as gene affiliation, with the NCBI2R package [30] which annotates lists of SNPs with current information from NCBI. We considered only SNPs located within a single gene in order to form gene groups without overlap. We focused our analysis on a list of 29 genes previously identified as having a main effect in GWAS.

The three methods tested yield different results, and only the PLS and G-GEE methods identify interactions. PCA detects only the main effect HLA-B and identifies no interactions. PLS detects the main effect HLA-B, but also identifies one interaction effect between the genes EOMES and BACH2. Our method G-GEE does not detect any main effects, but it shows two significant interactions, the first between the genes HLA-B and SULT1A1 and the second between IL23R and ERAP2.

Thyroid carcinomas

Thyroid cancers are thought to be related to a number of environmental and genetic predisposing factors and can be classified in various types and subtypes. Most association studies have focused on main effects but only a limited number of genes were identified. Recently, some papers focus on the detection of epistatic interactions [23, 31]. We applied our proposed approach on the two data sets used in Luzón-Toro et al. [23] regarding two rare tumours, sporadic medullary thyroid carcinoma (sMTC) and juvenile papillary thyroid carcinoma (jPTC). Affymetrix Genome-Wide Human SNP 6.0 arrays were used to hybridized DNA. The data set related to sMTC contains 66 cases and the jPTC data set 30 cases. The same 125 healthy controls and 232, 607 SNPs were used for both studies. As for the ankylosing spondylitis dataset, we obtained gene affiliation for each SNP with the NCBI2R package and considered only SNPs located within a single gene. We focused the analysis of the sMTC data set on a list of 10 genes, 3 of these genes (CHFR, AC016582.2 and C8orf37) were chosen following the conclusions of Luzón-Toro et al., the others because they contained markers that were susceptible to be associated with the disease from univariate analysis. The analysis of the jPTC data set was realized on a list of 20 genes among them we can cite DIO, RP11-648K4.2, LOXL1, DMGDH, PAX8 and STK17B from which epistatic interaction were already detected (even if the interaction between PAX8 and STK17B was identified in a study concerning papillary thyroid and not the juvenile form). The 14 others genes contained susceptible associated individual markers from univariate analysis. Regarding the sMTC study, G-GEE identifies one interaction between genes NCK1 and TRIQK. PCA detects only one main effect for the gene TRIQK whereas none effects were identified with PLS. Concerning the jPTC data set, 3 interactions were identified by G-GEE (NCAM1 and MNDA, MNDA and STK17B, LOC105370481 and STX3). PLS identifies 2 interactions (LOC105370236 and LOC105370481, LOC105370236 and PIKFYVE) and PCA detects only one main effect for the gene LOC105370481. We note that the effects detected with our approach concerned different genes that the ones identified in the presented previous studies (except for the gene STK17B). More analysese are needed to better understand these differences.

Inflammatory bowel disease

Although the etiology of Inflammatory Bowel Disease (IBD) is not completely understood, previous studies have underlined the contribution of an important genetic susceptibility. Recently, Martinez-Chamorro et al. [32] detected an epistatic interaction between the genes NOD2 and TLR4. We applied our approach to the Wellcome Trust Case-Control Consortium genome-wide association dataset for Inflammatory Bowel Disease. The data contains 1949 case for 159 960 SNPs genotyped by Affymetrix. The control group was constituted of 1972 individuals from the Wellcome Trust Case-Control Consortium genome-wide association dataset for hypertension. As for the two previous real data analysis, we obtained gene affiliation for each SNP with the NCBI2R package and considered only SNPs located within a single gene. The analysis was realized on a list of 22 genes that contain SNPs that are suspected to be associated with IBD from an univariate analysis. The two genes NOD2 and TLR4 were added to the list as they were previously detected as having an epistatic interaction. G-GEE identifies one interaction between the genes LOC105376008 and CACNB2 whereas PCA detects 9 main effects (IL23R, PODN, ATG16L1, C5orf56, DNAH11, LOC105378282, HSD17B12, LINC00558, ADCY4) but none interaction. Finally PLS identifies 3 main effects for the genes IL23R, PODN and DNAH11 as well as 2 interactions the first one between the genes PODN and FCRLA, the second one between PVT1 and NOD2.

Discussion

The results obtained in both simulation studies point to a certain confusion between main and interaction effects. When simulated interaction and main effects involve different genes, the methods tend to detect as interaction effects the pairs of genes simulated to have main effects and, conversely, to detect as main effects the genes simulated to having interaction effects.

Overall, G-GEE tends to detect more false interactions than false main effect whereas PLS and PCA tend to detect more false main effects though PLS tends to attribute a false interaction effect between genes 1 and 2. This type of confusion may explain the U-shaped power curve for G-GEE observed in the first simulation study (Fig. 1 a). As the problem becomes harder, the genetic effects of both genes are preferentially assigned to the interaction effect, implying a better power to detect interaction where R ² values are small. Finally, we remark that for G-GEE false detections of main effects are more frequent when the PCA phenotype simulation model is used, whereas for the PLS and PCA methods, where the number of false detections for main effects is higher when the Wang Pathway phenotype simulation model is used.

Other observations regarding the power of the different methods can be made with these simulation results. PLS has more trouble than PCA and G-GEE in detecting interaction effects, and has a tendency to detect the first main effect with a higher power than the second main effect when two main effects are simulated. For all methods, the power to detect interactions increases more slowly with respect to R ² when simulations are performed using real data genotypes than with fully simulated genotypes, but we observe that in the first setting the curve representing the interaction power of G-GEE is detached from the others, reflecting the superior performance of G-GEE over PLS and PCA. Note that the power of G-GEE to detect main effects is always less than that of PCA and PLS when R ²<0.4 [see Additional file 1 and 5]. In short, G-GEE performs better when detecting interactions than when detecting main effects.

Conclusions

In this paper we compared different approaches for modelling gene-gene epistasis in a penalized regression framework. Our primary concern was the detection of interaction effects, and for this purpose we defined a general model and tested different interaction terms. We focused our analysis at the gene scale and compared three ways to design the interaction term. Some methods were inspired by previous proposed approaches based on dimensional reduction methods including Principal Component Analysis (PCA) and Partial Least Square analysis (PLS). We additionally proposed a new interaction modeling approach that we called Gene-Gene Eigen-Epistasis (G-GEE), where one interaction variable is built for each couple. The interaction variable was defined based on a criterion that maximizes the covariance between the phenotype and the pairwise SNP product matrix of the two genes. The interaction components were then introduced in a Group Lasso penalized regression model that takes the gene structure into account and is capable of handling a large number of genes simultaneously.

A power study of the different methods based on two different simulation schemes (simplified and realistic) provided us with a rich body of information. Across various papers in the literature we find comparisons of similar methods that use different phenotype simulation settings. In the present work we compared two simulation models. Our first model was from a previous study [17] that simulated the interaction component of each couple in an SNP pairwise product fashion. Our second model defined the interaction component as a pairwise product of representative variables of each gene. Overall the G-GEE method performed well in detecting interactions in all the settings that were tested, although it was not always able to do so in the settings where main and interaction effects involved different genes. The power of the PCA method is highly dependent on the phenotype simulation model, because of the similarity between the second phenotype simulation model and the estimation model of the PCA method. The PLS method is characterized by a lack of power in detecting interactions. PLS performs well only when the related main effects are also present. When the simulated main and interaction effects do not concern the same genes, the detection capability of the PLS approach collapses dramatically.

For all methods we observed a confusion phenomenon when active genes are not simulated with both main and interaction effects. False detections of interactions concern genes that were simulated to have main effects, and false detections of main effects concern genes simulated to have interaction effects. This phenomenon reveals the difficulty that all methods encounter in clearly distinguishing the different types of effects. There are more false main detections when using methods such as PCA and PLS that are better at detecting main effects (except when the phenotype is simulated using the PCA model). As for interaction effects, the G-GEE methods make more false interaction detections than PCA and PLS.

When genotypes are fully simulated in the simplified simulation study, the G-GEE and PCA approaches performed better when the PCA phenotype simulation model was used, whereas the PLS method was not very sensitive to the choice of phenotype simulation model. Unlike PCA and PLS, G-GEE is better at detecting interaction effects than at detecting main effects when simulations use a real data set. Since the simulation study using realistic data is meant to mimic real genotype data structure, we conclude that in a real context G-GEE will be better at detecting interaction effects than main effects.

In comparison to SNP-SNP interaction approaches, the gene-scale dimension of our proposed method means that considerably fewer interaction variables need to be considered within a genetic region. This reduction in problem size allows larger problems to be handled. Moreover, a penalized regression method allows a true multivariate approach over a larger number of genes. It also extends other proposed gene-scale approaches, such as that presented by Wang et al. [8]. The ability to handle a relatively large number of genes simultaneously makes the detection of interactions between different genetic regions possible. This might be useful as an initial step, prior to using SNP-SNP interaction methods that may provide more accurate information.

As the G-GEE method is not able yet to consider all human genes at the same time, it is necessary to specify a list of genes to be explored for potential interactions. Given that its power to detect main effects is low, for the detection of main effects it will be safer to use previously acquired knowledge of the genetic effects, or to use a pre-processing method. Another limitation of the method is gene size. Computing the gene Eigen-Epistasis vector for two genes of size p _r and p _s requires an n×(p _r p _s ) matrix to be computed.

Prospects are improving the G-GEE method’s performance by optimizing the computational cost and exploring new interaction functions to be plugged into the G-GEE criterion.