1 Introduction

Complex human diseases are determined not only by genetic variants, but may also be affected by environmental factors and the interplay between them. Changes in gene expression under different environmental conditions reveal the interaction between genes and the environment. These changes are less likely due to changes in the gene sequence itself, but to structural changes such as DNA methylation or histone modification that consequently play a regulatory role and modulate gene expression. Such epigenetic changes have been increasing recognized as the epigenetic basis of gene-environment (G×E) interaction (Liu, Li & Tollefsbol, 2008). Identification of G×E interaction could shed novel insights into the phenotypic plasticity of complex disease phenotypes (Feinberg, 2004).

In a typical G×E interaction study, the environmental factor can be either discrete or continuous. For example, smoking can be a discrete variable when evaluating the risk of asthma. When environmental variables are measured on a continuous scale, a clearer picture of the interaction can be assessed since the varying patterns of genetic effects responsive to environmental changes can be traced, leading to a better understanding of the genetic heterogeneity under different environmental stimuli (Ma et al., 2011; Wu & Cui, 2013). As illustrated in Wu and Cui (2013), one can assess the nonlinear G×E interaction when an environmental factor is measured in a continuous scale. For example, individual obesity can be a factor when evaluating the risk of hypertension. One can assess the nonlinear effect of a genetic factor on the risk of hypertension considering the heterogeneity of individual obese conditions in a population, leading to a better understanding of disease heterogeneity.

When assessing G×E interactions, investigators have focused predominantly on single variant based analysis, such as the parametric methods in Guo (2000), the semi-parametric methods in Chatterjee and Carroll (2005), Chen, Chatterjee, and Carroll (2013), and Maity et al. (2009), and the non-parametric methods in Ma et al. (2011) and Wu and Cui (2013). Recently, there has been a significant increase in set-based genetic association studies focusing on a set of variants, for example, the gene-centric analysis of Cui et al. (2008), the gene-set analysis of Schaid et al. (2012) and Efron and Tibshirani (2007), and the pathway-based analysis of Wang, Li, and Hakonarson (2011). By assessing the joint function of multiple variants in a set, one can obtain a better interpretation of the disease signals and gain novel insights into disease etiology. Motivated by these set-based association studies, we propose a set-based framework to investigate how variants in a gene-set moderated by an environment factor affect disease and in what form.

In a typical set-based association study, the number of variants d within a genetic system can be relatively large compared to the sample size, which makes the regression coefficients estimation instable. The problem can be approached from the perspective of variable selection. In this work, we extend our previous work on nonlinear gene-environment interaction from a single variant based analysis to a multiple variant based analysis under a penalized regression framework. We include variants that belong to a particular gene-set or pathway which potentially interact with one or multiple environment factors through an additive varying-coefficient model. We propose to select genetic variants with coefficient functions that are varying, non-zero constant and zero corresponding to cases with G×E interactions, no G×E interactions and no genetic effects, respectively. Our approach employs the power and merits of variable selection by simultaneously fitting all the variants in a genetic system into a regression model, therefore avoiding the limitation of multiple testing corrections, especially when the data dimension is large.

This paper is organized as follows. In Section 2, we describe the penalized least square estimation procedure via B-spline basis expansion and smoothly clipped absolute deviation (SCAD) penalty, as well as the computational algorithms. We also present the theoretical results including consistency in variable selection and show the optimal convergence rates of the estimates of varying effects. We show that the estimates of non-zero constant coefficients enjoy the oracle property in the sense that the asymptotic distribution of the non-zero constant coefficient function is the same as that when the true model is known a priori. The merit of the proposed method is demonstrated through extensive simulation studies in Section 3 and real data analysis in Section 4. The technical proofs are relegated to the A.

2 Methods

3 Simulation

The performance of the proposed method was demonstrated through extensive simulation studies. We used the percentage of choosing the true model out of total R replicates, defined as the oracle percentage, to evaluate the accuracy of variable selection by identifying varying, non-zero constant and zero effects. The precision of estimation was assessed by integrated mean squared error (IMSE). Let β^j(r) be the estimator of a nonparametric function β_j in the rth (1 ⩽ r ⩽ R) replication, and {zm}m=1ngrid be the grid points where β^j(r) was evaluated. We used the integrated mean squared error (IMSE) of β^k(z), defined as

to evaluate the estimation accuracy of coefficient β_j, and the total integrated mean squared error (TIMSE) of all the d coefficients, defined as TIMSE=∑j=1dβ^j(z), to evaluate the overall estimation accuracy. Note that IMSE(β^_j) is reduced to MSE(β^_j) when β^_j is a constant. The percentage of correctly selecting each individual true functions (defined as the selection ratio) was used to evaluate the selection performance.

We considered multiple genetic factors X obtained from a gene-set or pathway, with the following additive VC model,

where SNP X_i’s were coded with 3 categories (1, 0, −1) for genotypes (AA, Aa, aa) respectively. We simulated the SNP genotype data based on the pairwise linkage disequilibrium(LD) structure. Suppose the two risk alleles A and B of two adjacent SNPs have the minor allele frequencies (MAFs) p_A and p_B, respectively, with LD denoted as δ. Then the frequencies of four haplotypes can be expressed as pab=(1−pA)(1−pB)+δ, pAb=pA(1−pB)−δ, paB=(1−pA)pB−δ, and pAB=pApB+δ. Assuming Hardy-Weinberg equilibrium, the SNP genotype at locus 1 can be simulated assuming a multinomial distribution with frequencies pA2, 2pA(1−pA) and (1−pA)2 for genotypes AA, Aa, aa, respectively. We can then simulate genotype for locus 2 based on the conditional probability. For example, P(BB|AA)=pAB2/pAA, P(Bb|AA)=pABpAb/pAA and P(bb|AA)=pab2/pAA. So conditional on genotype AA at locus 1, the genotype at locus 2 can be generated according to a multinomial distribution with the derived probabilities. The advantage of this simulation is that we can control the pairwise LD structure between adjacent SNPs. We assumed pairwise correlation of r = 0.5 which leads to δ=r(pA(1−pA)pB(1−pB)). To save space, we omitted the detailed simulation information which can be found in Cui et al. (2008). The coefficient functions were set as: β1(z)=sin⁡(2πz), β2(z)=2−3cos⁡{(6z−5)π/3}, β3(z)=3(2z−1)3, β4(z)=2, β5(z)=2.5, and βj(z)=0 for j > 5. We evaluated the performance under n = 500 with 500 replicates. Better performance results for large samples (n > 500) were observed, but were omitted to save space.

Figure 1 shows the selection ratio when d = 10, under different combinations of MAF and error distribution. The height of the bars represents the selection percentage out of 500 replicates. The selection performance is better under the normal error distribution, with relatively higher selection rate for the first five true functions and lower false selection ratio for the rest, compared to the results obtained under the t(3) error. In genetic association studies, model performance generally improves as the MAF increases. The same trend is observed under our variable selection framework. For example, a higher false selection ratio was observed under the t(3) error when p = 0.1. The false selection ratio decreases as MAF increases to 0.3. The result for d = 50 is presented in Figure 2, which shows a very similar pattern. The results demonstrate the stable performance of the proposed variable selection method.

Figure 1:

The selection ratio under different error distributions for different coefficient functions when d = 10. The horizontal axis represents the SNPs.

Table 1 lists the oracle percentage (%) of choosing the true model out of all the simulation replicates, the IMSE (inside the panel), and TIMSE (the last row) for the case with d = 10. In general, the model selection performance improves as the MAF increases from 0.1 to 0.5. For example, the oracle percentage increases from 0.72 to 0.91 under the t(3) error with SCAD penalty, when the MAF increases from 0.1 to 0.3. We observed dramatic reduction on the IMSE and TIMSE as the MAF increases. Under the normal error, the TIMSE is 0.4205 which reduces to 0.2007 when the MAF increases to 0.3 and further reduces to 0.1895 when p = 0.5. This result is consistent with the general observation in a genetic association study in which typically a model performs better as the MAF increases. It is worth mentioning that we observed dramatic improvement in model performance when the MAF increases from 0.1 to 0.3, compared to the improvement when the MAF increases from 0.3 to 0.5. For example, the IMSE for β2(u) reduces from 0.3285 to 0.1600, a 51% reduction when p increases from 0.1 to 0.3, while there is only a 1% reduction when p increases from 0.3 to 0.5 under the t(3) error distribution for the SCAD penalty. This empirical observation shows the stable performance of the model under moderate allele frequency.

Another observation from the simulation is that the model performs better under the normal error than under the t(3) error. We observed a larger oracle percentage, smaller IMSE and TIMSE for the coefficient functions under the normal error compared to the t(3) error. For example, the TIMSE for the SCAD penalty is 0.4205 under the normal error, while it is 2.9342 under the t(3) error for fixed p = 0.1. In addition, the oracle percentage, IMSE and TIMSE under the normal error are all quite similar as those obtained as if the truth were known (the oracle) in all cases, demonstrating the stable selection performance of the SCAD penalty.

A similar pattern was observed when the data dimension increases to 50 (Table 2). As the MAF increases from 0.1 to 0.3, we observed sharply decreased IMSE and TIMSE. Compared to the low dimensional case when d = 10, the performance under p = 0.1 is relatively unstable. For example, the TIMSE for the SCAD method is 3.3644 when d = 50, compared to 0.4205 when d = 10 under the normal error and p = 0.1. However, we observed dramatic reduction in TIMSE when the MAF increases to 0.3 under d = 50. Thus, one has to be very careful about the interpretation of the selection result under low MAF in real data analysis. We did additional simulations when the sample size increases to 1000 and observed consistently improved results under different scenarios (data not shown). In summary, the SCAD penalty function shows consistently good performance and can separate varying, constant and zero effects under moderate allele frequencies. Coupled with the results shown in Figure 1 and Figure 2, the proposed variable selection method shows relatively stable performance to assess gene-environment interactions.

Figure 2:

The selection ratio under different error distributions for different coefficient functions when d = 50. The horizontal axis represents the SNPs.

To further assess the false positive controls of the proposed method, we generated the response from the intercept only model, i.e. Yi=β0(Zi)+εi. There are no main and interaction effects associated with the disease phenotype. The average number of false positive effects for (d,error)=(10,N(0,1)),(10,t(3)),(50,N(0,1)),(50,t(3))) setups are 0.004, 0.042, 0.002 and 0.036, respectively. Overall, the proposed method achieves satisfactory false positive controls under the null model.

4 Case study

The body mass index of the mother (MBMI) is often used as a measure of the mothers’ body shape and degree of obesity. Since the baby resides inside its mother’s womb, its environment is defined through its mother. Increasing evidence indicates that both pre-pregnant weight (BMI) and weight gain in pregnancy have a major influence on babies’ birth weight (Stamnes Koepp et al., 2012). Due to the complicated interaction between the genes of the fetus and the mother’s level of obesity, the birth weight might be different for a fetus with the same genes but under different environment conditions. Thus, variation in birth weight could be totally or partially explained by the underlying genetic machinery and how those genes respond to the mother’s obesity to affect birth weight.

We applied the method to a real dataset from a study conducted in the Department of Obstetrics and Gynecology at Sotero del Rio Hospital in Puente Alto, Chile. The initial objective of the study was to pinpoint genetic variants associated with a binary response indicating large for gestational age (LGA) or small for gestational age (SGA) infants based on the birthweight of new born babies. After data cleaning by removing SNPs with MAF less than 0.05 or deviation from the Hardy-Weinberg equilibrium, the dataset contains 1536 new born babies genotyped with 189 candidate genes covering 660 single nucleotide polymorphisms (SNPs).

Genes were mapped to the KEGG pathway using the GATHER software which can be accessed at http://gather.genome.duke.edu. A total 30 pathways based on 189 candidate genes were retrieved. We treated the mother’s BMI as the environmental factor and the baby’s birth weight as the response variable; this was standardized before fitting to the model. Since some genes were mapped to multiple pathways, we did the variable selection for each pathway separately. Table 3 shows the selection results with SNP ID, the gene and pathway name the SNP(s) belong(s) to and the selected effect. Two SNPs in gene IL2 were mapped to two pathways and both SNPs consistently show varying effects in the two pathways. SNP rs2069762 in gene IL2 was previously reported to be associated with preterm birth and low birthweight in a Japanese population study (Sata F et al., 2009). Several other SNPs in gene IL1B were also reported to be associated with low birthweight in that paper. In addition, one SNP in gene IL1B in the Toll-like receptor signalling pathway was selected as a varying effect. Two SNPs in gene COL1A2 were mapped to two pathways and both were selected as varying effects. SNP rs997049 in gene IL1R1 was selected as a constant effect in two different pathways.

Figure 3 plots the varying coefficient function for two SNPs, SNP rs2039762 in the Cytokine-cytokine receptor interaction pathway and SNP rs2301643 in the ECM-receptor interaction pathway. The varying pattern of the function over mother’s BMI indicates the nonlinear interaction of the SNPs with mother’s BMI condition to affect birth weight. When fitting a linear interaction model, no SNPs show significant interaction with mother’s BMI (data not shown).

Figure 3:

The estimated varying coefficient function for SNP rs2039762 in the Cytokine-cytokine receptor interaction pathway and SNP rs2301643 in the ECM-receptor interaction pathway.

5 Discussion

The significance of G×E interactions in complex human disease traits has stimulated widespread discussion. As reviewed in Cornelis et al. (2011), a number of statistical models have been proposed to assess gene effect under different environmental exposures. The success of gene set based association analysis, as shown in Wang, Li, and Hakonarson (2011), Cui et al. (2008), Wu and Cui (2013), and Schaid et al. (2012), motivated us to propose a high dimensional variable selection approach to understand the mechanism of G×E interactions associated with complex diseases. We adopted a penalized regression method within the VC model framework to investigate how multiple variants within a genetic system are moderated by environmental factors to influence the phenotypic response.

Within the model-based regression framework, most G×E interactions are modeled via a product term between a G and an E variable (Hutter et al., 2013), so the contribution of a genetic variant to the phenotypic variation is considered as a linear function in the environmental factor. Any non-linear interaction can be pursued to relax the linearity assumption (Ma et al., 2011; Wu & Cui, 2013). As pointed out by one reviewer, statistical interactions introduced by R.A. Fisher, are defined as deviation from a generalized linear model, which implicitly suggests a nonlinear relationship and is more general. To avoid confusion with the nonlinear G×E interaction presented in this work, we make it clear that our nonlinear G×E interaction refers to the effect of a genetic variant assessed as a nonlinear function of an environment variable.

In a G×E study, people are typically interested in assessing variants which are sensitive to environment changes and those that are not. We can determine if a particular genetic variant is sensitive to environmental stimuli by examining the status of the coefficient function. Varying-coefficients and constants can be separated through B-spline basis expansions under a penalized framework. The varying coefficients correspond to G×E effects and the constant effects correspond to no interaction effects. Through another penalty function, we can further shrink the constant effect into zero if the corresponding SNP has no genetic effect. We developed a two-stage iterative estimation procedure with double SCAD penalty functions. Although the two-stage strategy has been adopted for regularized quantile regression with adaptive LASSO in Tang et al. (2012), our work significantly differs in that we focused on the regularized least square regression and rigorously establish the asymptotic properties of the nonconvex double SCAD estimator under suitable regularity conditions. The potential of non-convex penalty functions in investigating G×E interactions is far from fully understood or explored. As a representative non-convex penalty function, the SCAD is adopted mainly due to its nice oracle properties as stepping stones for building statistically sound and practically useful models to accommodate more complex data structures. It is worth to mention that our method is fundamentally different from the work of Xue and Qu (2012) and Antoniadis, Gijbels, and Lambert-Lacroix (2014) in which the authors developed a variable selection framework under the additive VC model to distinguish zero vs varying coefficients. They did not distinguish non-zero constant vs varying coefficients which is one of the key objects in understanding the mechanisms of G×E interaction. Identification of the constant coefficients in the varying coefficient models is closely related to the estimation of linear part in additive models. This line of work includes Hu and Xia (2012) and Zhang, Cheng, and Liu (2011). None of the existing studies closely explore the automatic structure identification and separation of different effects under the G×E framework.

The current work only demonstrates the case with one environmental factor. It is broadly recognized that the etiology of many complex disease is less likely to be affected by one environmental factor but is more likely to be heterogeneous. When multiple continuously measured environmental factors (say K₁) are measured (denoted as Z₁), we can extend the current model to a more general case formulated as follows,

where X₀ = 1. The same estimation and variable selection framework can be applied to select important genetic players that show sensitivity to different environmental stimuli. When discrete environmental variables such as smoking status are also available, denote Z₂ as a collection of K₂ such variables, then we can fit the following model

the partial linear varying-coefficient model. In addition to the two penalty functions specified in this work, an additional penalty function should be imposed for {α}lj to select important variants showing interaction with Z₂. In case of a binary response, we are interested in modeling E[Y|Z1,Z2,X]=∑j=0d{∑k=1K1βkj(Z1k)+∑l=1K2αljZ2l}Xj. We will investigate this in future studies.

The proposed method is not only restricted to quantitative phenotypes and can be extended to other types of phenotypes. For example, in cancer prognostic studies, it can be modified as the Accelerated failure time (AFT) model to accommodate the survival outcomes. Binary phenotypes significantly differ from quantitative and survival outcomes in that they contain much less information, hence the accuracy of estimating nonlinear interactions might be sacrificed. Nevertheless, extension to the binary case can be done by developing a coordinate descent (CD) based iteratively reweighted least squares (IRLS) algorithm under the regularized logistic regression framework. The CD based IRLS algorithm have been extensively used to extend regularized variable selection methods from continuous phenotypes to binary phenotypes such as in case control studies.

In the model, we did not include any covariates. However, the proposed varying coefficient model can be readily modified to allow for covariate effects. Typically, the covariates included in the model are predetermined as important ones and are in low dimension, so their effects are not subject to penalization. Assuming there are no interactions between genes and those covariates (those with interactions will be included in the model), one can fit a regression model by regressing Y against those covariates only, assuming either linear on nonlinear effects, Then focusing the obtained residuals (after removing the covariates effects) to do the rest of the analysis by fitting the models described in this paper. It is also worth mentioning that the real data analysis in this work does not take other covariates (e.g. gender and mother’s gestational age) effects into account, which may lead to biased results. Due to this limitation, readers should be cautious when interpreting the real data analysis results.