An integrative analysis of genomic and exposomic data for complex traits and phenotypic prediction.

Complementary to the genome, the concept of exposome has been proposed to capture the totality of human environmental exposures. While there has been some recent progress on the construction of the exposome, few tools exist that can integrate the genome and exposome for complex trait analyses. Here we propose a linear mixed model approach to bridge this gap, which jointly models the random effects of the two omics layers on phenotypes of complex traits. We illustrate our approach using traits from the UK Biobank (e.g., BMI and height for N ~ 35,000) with a small fraction of the exposome that comprises 28 lifestyle factors. The joint model of the genome and exposome explains substantially more phenotypic variance and significantly improves phenotypic prediction accuracy, compared to the model based on the genome alone. The additional phenotypic variance captured by the exposome includes its additive effects as well as non-additive effects such as genome-exposome (gxe) and exposome-exposome (exe) interactions. For example, 19% of variation in BMI is explained by additive effects of the genome, while additional 7.2% by additive effects of the exposome, 1.9% by exe interactions and 4.5% by gxe interactions. Correspondingly, the prediction accuracy for BMI, computed using Pearson's correlation between the observed and predicted phenotypes, improves from 0.15 (based on the genome alone) to 0.35 (based on the genome and exposome). We also show, using established theories, that integrating genomic and exposomic data can be an effective way of attaining a clinically meaningful level of prediction accuracy for disease traits. In conclusion, the genomic and exposomic effects can contribute to phenotypic variation via their latent relationships, i.e. genome-exposome correlation, and gxe and exe interactions, and modelling these effects has a potential to improve phenotypic prediction accuracy and thus holds a great promise for future clinical practice.

Introduction

Both genetic and environmental factors underlie phenotypic variance of complex traits. Understanding the influences of these factors not only helps explain why individuals differ from one another in phenotypes but also helps predict future phenotypes, such as disease diagnoses. The proliferation of genotypic data in the past decades, along with developments in relevant analytic tools, have already contributed a great deal to understanding phenotypic variations of complex traits[1-9], and enabled phenotypic predictions at a level of accuracy for potential use in clinical settings[10-12]. However, these understandings and predictions are bounded by the heritability of the traits, and for many complex traits, large phenotypic variation remains unexplained, suggesting substantial environmental contributions to phenotypic variance.

Complementary to the genome, the concept of exposome has been proposed to capture the totality of human environmental exposures, encompassing external as well as internal environments over the lifetime of a given individual[13-15]. Similar to genotypes, exposomic variables are standardised across cohorts[16]. Since the inception of the concept, considerable efforts have been made to assess and characterise the exposome[17]. For example, the Human Early-Life Exposome project is a European collaborative effort established to characterize the early-life exposome which includes all environmental hazards that mothers and children are exposed to[18]. Despite the progress in the construction of the exposome, few analytic tools exist to date that can integrate genomic and exposomic data for complex trait analyses. We hypothesize that exposomic variables do not only affect phenotypes on their own but also interact among each other[19,20] and with genotypes[20,21]. In addition, the estimation of exposomic effects and genomic effects on phenotypes could be biased, if these effects are correlated but the estimation model assumes otherwise[22]. Hence, tools that integrate genomic and exposomic data are required to capture variance as well as covariance components of phenotypes.

Here we propose a versatile linear mixed model that fulfils these requirements. The proposed approach jointly models the random effects of the genome and exposome and can be extended to capture genome-exposome and exposome-exposome interactions and genome-exposome correlations in the phenotypic analysis of a complex trait. It also allows us to model exposomic effects modulated by one or a few specific environmental variables. We demonstrate the proposed approach using traits from the UK biobank with 11 complex traits and 28 lifestyle exposures that were measured using a standard protocol.

Results

Method overview

We used a novel linear mixed model (LMM) to jointly model the effects of the genome and exposome on the phenotypes of a complex trait. The exposome here is restricted to 28 lifestyle exposures that were measured using a standard protocol (see “Methods”). Our model has three key features. First, it allows estimation of the correlation between genomic and exposomic effects, relaxing the assumption of independence between those effects as in a conventional LMM[22]. Second, the model can capture both additive and non-additive effects of the exposome and genome, i.e. pairwise interactions between exposomic variables (exe interactions; e.g.[19]) and interactions between exposomic variables and genotypes (i.e., gxe interactions; e.g.[21]). Third, the model can handle correlated exposomic variables (see Simulations in “Methods”) that may cause biased variance estimations of exposomic variables (e.g.[20]).

To illustrate the use of the model with real data, we selected 11 complex traits from the UK Biobank with heritability estimates above 0.05, including BMI, sitting height and years of education etc. (https://nealelab.github.io/UKBB_ldsc/), along with 28 lifestyle variables, including alcohol use, smoking, physical activity and dietary composition (see “Methods” for a detailed description). We performed the following analyses. First, for each trait, we used various models to estimate variance components of the additive and non-additive effects of the exposome and genome, including exe interactions and gxe interactions. The significance of the variance components was determined through a series of model comparisons using likelihood ratio tests (Table 1). Second, we extended the proposed model to examine the extent to which exposomic effects are modulated by covariates such as age, sex and socio-economic status (i.e., exc interactions). Third, we used fivefold cross validation to show that the prediction accuracy increased significantly after accounting for the exposomic effects and exe interactions. Finally, we explored the potential clinical use of the proposed integrative analysis of genomic and exposomic data, by projecting its prediction accuracy for a disease trait in terms of the area under the receiver operating characteristic curve (AUC). The projection was based on well-established theories[23-30] that express AUC as a function of sample size, proportions of variance explained by genomic and exposomic effects and the population prevalence of the disease.

Table 1

P-values for estimated variance components of selected traits.

Significant estimates are orange-coloured after applying Bonferroni corrected alpha level for each model comparison = 0.05/66 = 7.6E−04.

See Table 2 for details of statistical models.

Table 2

Model equations and their assumed sample variance–covariance matrices.

	Model equation	Matrix notation	Sample variance–covariance matrix
	For individual i = 1, 2, …., n,	For y = (y1, y2, …., yn),
(i)	where yi is the phenotype, μ is the grand mean, aij is the standardized SNP genotype at locus j, m1 is the total number of SNPs, αj is the random effect of the SNP that is assumed to be normal with mean zero and variance σg2/m_1, and εi is the residual assumed to be normal with mean zero and variance σε2.	y=μ1_n+g+ε where g=Aα^t=a_11⋯a_1m_1⋮⋱⋮a_n1⋯a_nm_1α_1⋮α_m_1	Where A is a n x m1 matrix that contains column-standardized genotypes (see the Matrix Notation), and I is the n x n identity matrix
(ii)	where bk is the kth exposomic variable, m2 is the total number of exposomic variables, and βk is the random effect of the exposomic variable that is assumed to be normal with mean zero and variance σe2/m_2. To avoid estimation bias due to multicollinearity, bk is transformed using a principal component analysis (see “Methods”)	y=μ1_n+g+e+ε where e=Ωβ^t=b_11⋯b_1m_2⋮⋱⋮b_n1⋯b_nm_2β_1⋮β_m_2	Where Ω is a n x m2 column-orthogonal and column-standardised matrix that contains the transformed exposomic variables multiplied by their right singular vectors (see ‘Principal component-based transformed variables for E’ in “Methods”)
(iii)	y_i=μ+g_i+e_i+ε_i	y=μ1_n+g+e+ε	σg2G+σe2E+G+E^t+G+E^t^tσ_ge+σε2I where G and E^t are the Cholesky decompositions of G and E^t, respectively, and σ_ge is the covariance between g and e
(iv)	where c_q is the qth pairwise interaction term between SNP genotypes and exposomic variables, and γ_q is the effect of the qth interaction term.γ_q is assumed to be normally distributed with mean zero and variance σg×e2/Q, and Q is the total number of interaction terms (Q=m_1m_2)	y=μ1_n+g+e+g×e+ε where g×e=Cγ^t=c_11⋯c_1Q⋮⋱⋮c_n1⋯c_nQγ_1⋮γ_Q, and C can be derived using the following pseudo-code with A=a_1⋯a_m_1; Ω=b_1⋯b_m_2; C=c_1⋯c_Q, and q = 1, 2 … Q for i = 1 to m1 { for j = 1 to m2 { c_q=a_i⊙b_j} }	σg2G+σe2E+σg×e2Γ+σε2I where Γ is a n x n matrix derived by the Hadamard product of G and E ,i.e., Γ=G⊙E=CC^t/(m1∗m2)
(v)	where x_p is the pth pairwise interaction term between exposomic variables, and when the two exposomic variables are identical, the interaction term becomes the quadratic term of the exposomic variable; θ_p is the effect of the pth interaction term and is assumed to be normally distributed with mean zero and variance σe×e2/P, and P is the total number of interaction terms (P = m2 (m2 + 1)/2). To avoid estimation bias due to multicollinearity, xp is transformed using a principal component analysis (see “Methods”)	y=μ1_n+g+e+e×e+ε where e×e=Xθ^t=x_11⋯x_1P⋮⋱⋮x_n1⋯x_nPθ_1⋮θ_P, and X can be derived using the following pseudo-code with Ω=b_1⋯b_m_2; X=x_1⋯x_P, and p = 1, 2 … P for i = 1 to m2 { for j = i to m2 { x_p=b_i⊙b_j} }
(vi)	y_i=μ+g_i+e_i+g×e_i+e×e_i+ε_i	y=μ1_n+g+e+g×e+e×e+ε	σg2G+σe2E+σg×e2Γ+σe×e2Θ+σε2I
(vii)	where λ_kl is the random effect of kth exposomic variable, b_k, modulated by the lth covariate c_l. λ_kl is assumed to be normally distributed with mean zero and variance σe_l2/m_2	y=μ1_n+g+e+∑l=1Le×c_l+ε where e×c_l is a n × 1 vector that can be derived by e_l⊙c_l, and e_l=b_11⋯b_1m_2⋮⋱⋮b_n1⋯b_nm_2λ_1l⋮λ_m_2l	σg2G+E⊙ϕKϕ^t+σε2I where ϕ=1_nc_1c_2⋯c_L and K=σe_02⋯σ_e_0e_L⋮⋱⋮σ_e_0e_L⋯σe_L2

Exposomic effects on phenotypes

In line with previous estimation (https://nealelab.github.io/UKBB_ldsc/), we found significant SNP-based heritability for all selected traits, with estimates ranging between 0.08 (years of education) and 0.52 (standing height; Fig. 1). We detected significant additive effects of the lifestyle-exposome on phenotypes of all traits (see Fig. 1 for e and Table 1 for p-values under H0 ). The magnitude of these additive effects, however, varied across traits. For example, the exposome accounted for 8.5% of the phenotypic variance of waist circumference, but less than 2.5% for height, standing height, heel bone mineral density and fluid intelligence. Importantly, the additive exposomic effects were mostly uncorrelated with the genetic effects (see Table 1 for p-values under H0 ; see Supplementary Table 1 for covariance estimates), which was notably different from the genome-transcriptome correlation[22].

Figure 1

Breakdown of phenotypic variance by the model with the best fit. The best model for each trait is derived from model comparisons shown in Table 1. g: additive genetic effects on phenotypes, e: additive effects of exposomic variables, exe: interaction effects between exposomic variables, gxe: interaction effects between genotypes and exposomic variables. Variance components are expressed as percentage of total phenotypic variance.

The estimated variance component of non-additive effects of the lifestyle-exposome (exe) was highly significant for 7 out 11 traits (Table 1), although they only account for ~ 1% to 2% of phenotypic variance (See Fig. 1 & Supplementary Table 2). By contrast, significant gxe interactions are only evident for BMI, weight and years of education (Table 1), but they could account for up to 9% of total phenotypic variance (years of education; Fig. 1 and Supplementary Table 2). The low presence of gxe signals is probably due to relatively low power of detecting gxe interactions, which is caused by a large number of pairs of gxe interaction terms to be estimated in the model, i.e. 28 (number of exposomic variables) × 1.3 million (number of SNPs) in this study. In addition, the identified gxe and exe interactions are largely independent to each other. This is evidenced by that both gxe and exe remained significant when being jointly modelled (see p-values under H0 and under H0 ).

By extending the proposed model to a reaction norm model (RNM; see “Methods”), we examined whether the additive exposomic effects on phenotype vary depending on specific covariates, which would be evidenced by the presence of significant exc interactions. Using single-covariate RNMs, we identified several significant exc interactions (Supplementary Table 3), noting that most covariates are lifestyle related, which are in line with the exe interactions found above. For each trait, we then fitted an RNM model that simultaneously includes all significant exc interactions identified from single-covariate RNM analyses. The variance estimates of exc interactions from the joint analyses are presented in Supplementary Table 4.

It is important to note that the estimation of exposomic effects is sensitive to the correlation structure of exposomic variables. Specifically, multicollinearity between exposomic variables would bias the estimate of (see simulations in “Methods”); and by extension, correlated exe interaction terms and gxe interaction terms (model equations iv and v in Table 2) could bias the estimates of and , as empirically observed in the simulations (see “Methods”). Without knowing the true values of variance components, transforming exposomic variables and interaction terms using a principal component analysis (see “Methods”) seems necessary prior to model fitting in order to avoid estimation bias due to multicollinearity. While transforming the exposomic variables and the exe interaction terms are computationally trivial, transforming the gxe interaction terms is computationally infeasible (28 × 1.3 million variables). Nonetheless, the variance of gxe interactions is small in general, suggesting that using the gxe interaction terms without the transformation (i.e., derived from in model equation iv of Table 2) is generally free from the estimation bias due to multicollinearity. Note that the largest variance estimate of gxe interactions in this study is ~ 0.09.

where yi is the phenotype, μ is the grand mean, aij is the standardized SNP genotype at locus j, m1 is the total number of SNPs, αj is the random effect of the SNP that is assumed to be normal with mean zero and variance , and εi is the residual assumed to be normal with mean zero and variance

where

Where is a n x m1 matrix that contains column-standardized genotypes (see the Matrix Notation), and is the n x n identity matrix

where bk is the kth exposomic variable, m2 is the total number of exposomic variables, and βk is the random effect of the exposomic variable that is assumed to be normal with mean zero and variance To avoid estimation bias due to multicollinearity, bk is transformed using a principal component analysis (see “Methods”)

where

Where is a n x m2 column-orthogonal and column-standardised matrix that contains the transformed exposomic variables multiplied by their right singular vectors (see ‘Principal component-based transformed variables for E’ in “Methods”)

where and are the Cholesky decompositions of and , respectively, and is the covariance between and

where is the qth pairwise interaction term between SNP genotypes and exposomic variables, and is the effect of the qth interaction term. is assumed to be normally distributed with mean zero and variance , and is the total number of interaction terms ()

where and can be derived using the following pseudo-code

with ; ; , and q = 1, 2 … Q

} }

where is a n x n matrix derived by the Hadamard product of and ,i.e.,

where is the pth pairwise interaction term between exposomic variables, and when the two exposomic variables are identical, the interaction term becomes the quadratic term of the exposomic variable; is the effect of the pth interaction term and is assumed to be normally distributed with mean zero and variance , and P is the total number of interaction terms (P = m2 (m2 + 1)/2). To avoid estimation bias due to multicollinearity, xp is transformed using a principal component analysis (see “Methods”)

where ,

and can be derived using the following pseudo-code with ; , and p = 1, 2 … P

} }

where is the random effect of kth exposomic variable, , modulated by the lth covariate . is assumed to be normally distributed with mean zero and variance

where is a n × 1 vector that can be derived by

, and

where and

Validation of exposomic effects

Figure 2a shows the phenotypic prediction accuracy based on genetic data alone. Using fivefold cross-validation, we found that including additive (e) and non-additive effects (exe) of the exposome, which were significant in the discovery dataset, could improve the phenotypic prediction accuracy in the target dataset. In general, the larger the variance estimates, the greater the prediction improvements (Fig. 2b,c), which indicates that the additive effects of the exposomic variables and exe interactions are genuine. Similarly, we also validated the exposomic effects modulated by specific covariates, by showing that the larger the total variance estimates of exc interactions, the greater the improvement of predication accuracy (Fig. 3). The validated exc interactions would in part explain the phenotypic variance due to residual x covariate interactions found in our previous studies[31-33].

Figure 2

Exposomic variables contribute to phenotypic variance and improve phenotypic prediction accuracy. The prediction accuracy of a given model was computed using the Pearson’s correlation coefficient between the observed and the predicted by the model. For all panels, the least squares line with 95% confidence band is based on a linear model that regressed either prediction accuracies (a) or predication accuracy improvements (b,c) by a model on variance component estimates of the model. The p-value is for the t-test statistic (df = 7) under the null hypothesis that the slope of the regression line is zero. = phenotypic variance explained by additive effects of the genome; = phenotypic variance explained by additive effects of the exposome; = phenotypic variance explained by exposome-by-exposome interactions; and = total phenotypic variance. (a) Phenotypic prediction accuracies by the baseline model that uses genomic data alone, i.e., y = g + ε, where g = phenotypic effects of the genome and ε = residuals. The larger the genetic variance, the greater the prediction accuracy. (b) Additive effects of the exposomic variables (i.e., e) contribute to phenotypic variance and improve phenotypic prediction accuracy. The greater the additive effects, the larger the gain in phenotypic prediction accuracy. A prediction accuracy improvement (on the y-axis) was derived by subtracting the prediction accuracy of the model y = g + ε from that of the model y = g + e + ε. (c) Exposome-by-exposome interactions (i.e., exe interactions) contribute to phenotypic variance and further improve phenotypic prediction accuracy. The greater the variance estimate of exe interactions, the larger the gain in phenotypic prediction accuracy. A prediction accuracy improvement (on the y-axis) was derived by subtracting the prediction accuracy of the model y = g + e + ε from that of the model y = g + e + exe + ε.

Figure 3

Positive relationship between phenotypic variance explained by exposome-by-covariate (exc) interaction effects and prediction accuracy improvement. Prediction accuracy improvement is computed by subtracting the prediction accuracy of the model y = g + e + ε from that of a model with multiple covariates (see model equation vii in Table 2) that are shown to interact with the exposome in univariate exc interaction analyses. The least squares line with 95% confidence band is based on a linear model that regressed prediction accuracy improvement on phenotypic variance explained by exc interactions. The p-value is for the t-test statistic (df = 7) under the null hypothesis that the slope of the regression line is zero. Significant covariates included for each trait can be found in Supplementary Table 3.

By contrast, although gxe interactions contribute to the phenotypic variance of BMI, weight and years of education (Table 1), the contribution did not lead to significant gains in phenotypic prediction accuracy (Supplementary Fig. 1). This was most likely due to a lack of power. i.e. the size of discovery samples was insufficient to accurately estimate an extremely large number of parameters, i.e., best linear unbiased prediction (BLUPs) of gxe interaction effects[23,27,28,34]. This is further verified using simulations (see Simulations in “Methods” and Supplementary Fig. 2).

Given the sample sizes of the discovery data sets (~ 28,000), the prediction accuracies of the model y = g + ε for the selected traits are only between 1/3 and 1/2 of the theoretical maximums (i.e., square root of heritability; Supplementary Fig. 3). They can improve, in theory, by increasing the sample size of discovery sets (Supplementary Fig. 3); or, as shown in the above, by accounting for the additive effects of the exposome and exe interactions (Fig. 2b,c). To examine prediction efficiency of the latter, we projected the observed prediction accuracies of the models y = g + e + ε and y = g + e + exe + ε onto the theoretical trajectory of prediction accuracies of the model y = g + ε as a function of the sample sizes of discovery datasets (Supplementary Fig. 3). As such, the use of exposomic information could improve phenotypic prediction accuracy to the same extent as a 1.2 to 14-fold increase in sample size, depending on the significance of the exposomic effects and their interactions (Fig. 4). Given the substantial costs and efforts associated with increasing sample size, the improved predictive accuracy by the models y = g + e + ε and y = g + e + exe + ε are considerable, despite the fact that the proportion of phenotypic variance explained by the exposome is small (see the x-axis of Fig. 2b,c).

Figure 4

Additional sample size required for the model y = g + ε to achieve the same level of prediction accuracy as y = g + e + ε (blue) and y = g + e + exe + ε (red). nt = sample size of training (or discovery) datasets.

Quantification of clinical relevance

We quantified the clinical relevance of the proposed model by exploring its prediction accuracy for quantitative traits and disease traits. For quantitative traits, we expressed the prediction accuracy of the model y = g + e + ε (i.e., correlation coefficient between the true and predicted phenotypes) as a function of the sample size of the discovery dataset, variances explained by the genome and exposome, and effective numbers of (independent) SNPs and exposomic variables (see “Methods”), using previous theoretical derivations[27-30,34]. Based on the derived expression [Eq. (6)], we computed the expected prediction accuracies for the quantitative traits used in this study and found that they agreed well with the observed prediction accuracies from the fivefold cross validation (Supplementary Fig. 4). We then extended the derived expression to disease traits in terms of the area under the operative characteristic curve [AUC; see Eq. (10) in “Methods” for details] using well-established theories[23-26]. AUC is a gold-standard measure used to evaluate how well a prediction model discriminates diseased from non-diseased individuals. An AUC between 0.7 and 0.8 is considered acceptable, 0.8 to 0.9 excellent, and above 0.9 outstanding[35]. Figure 5 shows the expected AUC of the proposed integrative analysis of genomic and exposomic data for disease traits of different values of population prevalence (k), assuming different amounts of variance (on the liability scale) explained by the genome and exposome and discovery sample sizes. For simplicity, we use to denote the total variance in disease liability explained by additive effects of the exposome and exe interactions as a whole.

Figure 5

Expected prediction accuracy of the proposed integrative analysis of genetic and exposomic data for disease traits of different prevalence (k) and heritability (h2) at varying levels of total variance explained by the exposome () and sample size of the discovery dataset (N). Diseases are assumed to have a liability of mean zero and variance 1, and both h2 and are on the disease liability scale. Prediction accuracy is measured using the area under the receiver operating characteristic (ROC) curve, with 0.7 to 0.8 generally being considered acceptable, 0.8 to 0.9 excellent, and above 0.9 outstanding. The assumed effective number of chromosome segments and the number of exposomic variables are 50,000 and 28, respectively, which are based on the genomic and exposomic data used in this study. However, varying the number of exposomic variables from 28 to 100 does not have a notable effect on the expected area under the ROC curve.

When setting to 0—that is, using no exposomic information at all—varying the heritability of disease liability h2 from 0 to 0.3 improves AUC from 0.5 to ~ 0.6 when the sample size of the discovery set is 50 k. This is in contrast to a twofold improvement, from 0.5 to ~ 0.7, when the sample size is 500 k. Thus, genomic prediction accuracy heavily relies on sample size, such that for a disease trait with a moderate heritability, a clinically meaningful level of accuracy (AUC ≥ 0.7) may not be attainable unless the sample size of the discovery dataset is substantially large (≥ 500 k). On the other hand, the benefit of using exposomic information to disease prediction can be realised with a relatively small discovery sample. This is evidenced by that when setting h2 to 0 (i.e., using no genomic information at all), increasing the value of has the same effects on AUC whether using a discovery sample of 50 k or 500 k individuals. Importantly, AUC consistently improves with increasing in all scenarios (Fig. 5). Thus, incorporating exposomic data is not only an efficient but also an effective way of improving prediction accuracy based on genomic data alone. Taken together, genomic prediction accuracy for disease traits is constrained by sample size; with a relatively small sample at hand, a desired level of prediction accuracy may only be achieved by combining genomic and exposomic information.

Comparison with existing models

The key model parameters of the proposed integrative analysis of genomic and exposomic data (IGE) compared to existing linear mixed models that incorporate genetic and environmental effects on phenotypes are outlined in Table 3. In general, IGE offers thus far the most detailed partition of phenotypic variance.

Table 3

Comparisons of methods (software packages) on the genomic and exposomic analysis of complex traits.

IGE (proposed method): integrative analysis of genomic and exposomic data.

StructLMM & GxEMM are existing linear mixed models that incorporate genetic and exposomic effects on phenotypes.

●: the parameter is included in the model, but the parameter estimate is not provided by the software package.

●●: the parameter is included in the model, and the parameter estimate is provided by the software package.

v(g): additive genetic variance due to either a single SNP or all common SNPs (i.e., whole genome).

v(gxe): GxE variance due to either interactions of a single SNP or all common SNPs with multiple exposomic variables.

GWEIS: Genome-wide by environment interaction study. Using the SNP BLUP method, the software for IGE (mtg2 v2.18) provides allele substitution effects of SNPs across environments, their standard errors and p-values. The StructLMM software provides allele substitution effects and p-values for GxE interactions.

cov(g,e): covariance between genomic and exposomic effects on phenotypes.

v(e): variance due to additive effects of exposomic variables.

v(exe): variance due to exposome x exposome interactions.

v(exc): variance due to exposome x covariate (e.g., demographics) interactions.

bivariate or multivariate analyses: analyses that simultaneously involve two or more traits.

Both IGE and GxEMM[36] are whole-genome approaches to the estimation of heritability and gxe interactions, although IGE is considered more comprehensive and versatile, which models variances explained by additive effects of exposomic variables, by exposome × exposome interactions, and by exposome × covariate (such as demographics) interactions; and covariance between genetic effects and exposomic effects (Table 3). Further, bivariate or multivariate IGE (i.e., simultaneously including two or more traits) can be feasibly performed using mtg2 version 2.18 (https://sites.google.com/site/honglee0707/mtg2).

In contrast, StructLMM has been developed primarily for a genome-wide by environment interaction study (GWEIS)[20] that examines one SNP at a time with a focus on association tests (providing p-values) for G × E interactions between the SNP genotypes and multiple exposomic variables. Using the well-established SNP BLUP method[2,37,38], IGE can also provide GWEIS summary statistics, including estimated allele substitution effects of all SNPs across environments, their standard errors and p-values. Note that SNP BLUP implemented in IGE can model all SNP jointly (a whole-genome approach). Nonetheless, one of the main scopes of this study is to provide unbiased estimates of exposomic variances, e.g., that is common to both StructLMM and IGE (Table 3 and Supplmentary Note 1). Importantly, correlated exposomic variables would cause biased estimation of (Supplementary Table 5) unless they are transformed to independent variables via a principal component analysis (“Methods”). To our knowledge, this transformation has not yet been implemented in any existing methods including StructLMM. Using results from simulations, we show that estimates by StructLMM are prone to bias due to correlated environments (Supplementary Table 5). The other model parameters such as , , and cannot be estimated by StructLMM (Table 3).

Discussion

Using our approach, we demonstrate the importance of the exposome for understanding individual differences in phenotypes. Although the ‘exposome’ constructed in this study comprises only 28 lifestyle factors, when integrated with genomic data, it explained between 2 to 10% additional phenotypic variance and significantly improved phenotypic prediction accuracy to a level equivalent to a 1.2 to 14-fold increase in sample size. The additional phenotypic variance is not only from additive effects of the exposome but also from its non-additive effects (exe) and genome–exposome interactions (gxe). We expect that as the construction of the exposome continues to progress, more phenotypic variance will be explained and greater improvements in phenotypic prediction accuracy will be gained. This would be particularly promising for phenotypic analysis and prediction of traits with small to little heritability component, such as ovarian and colorectal cancer[39].

We noted that when exposomic variables are correlated, the variance estimate of additive effects of exposomic variables is biased unless these variables are transformed using a principal component analysis (i.e. E in Table 2 should be based on transformed variables). By extension, this would apply to exe interaction terms and gxe interactions terms, unless the proportions of phenotypic variance explained by these interaction effects are small (< 10%), as shown in our simulations. These observations have important implications for modelling environmental effects in LMMs. Recently, Moore et al.[20] proposed the structured linear mixed model (StructLMM) that incorporates random effects of multiple environments in order to study the interactions between these environments and genotypes of a single SNP (i.e., gxe interactions). However, the environmental variables in StructLMM are not transformed, even though they are very likely correlated, which would have biased the variance estimate of environmental effects. Consequently, it remains uncertain the extent to which the estimation bias affects the StructLMM-based test statistics for detecting gxe interactions.

Depending on the research question at hand, the construction of the exposome may be guided by causal analyses. A meaningful exposome may only contain causal information. Examples may include lifestyles that potentially alter the molecular pathways or the pathogenesis of the main trait, or biomarkers that potentially explain possible molecular pathways underlying the phenotypes. As a contrast, in our BMI analysis, for example, it is not useful to include weight and height as part of the exposome, even though they would explain a large amount of phenotypic variance. This is because variations in these traits inform nothing other than the fact that they are correlated with the trait.

Heritability estimates were slightly reduced after including more variance components (result not shown). We considered two possibilities. First, the exposome may mediate part of additive genetic effects on phenotypes. For example, some SNPs affect smoking status, which in turn affect BMI. A model that simultaneously includes genetic and exposomic data would account for smoking status and their genetic effects, and hence gives arise to reduced heritability estimates. Second, there is a genuine correlation between exposomic and genomic effects in some latent mechanism. It is noted that there are marginally significant correlation estimates, which were not significant after Bonferroni correction. Such correlation may be because people who have similar genotypes may somehow share similar exposures i.e. genotype-environment correlation[40].

In conclusion, the genomic and exposomic effects can contribute to phenotypic variation via their latent relationships, i.e. genome-exposome correlation, and gxe and exe interactions, for which our proposed method can provide reliable estimates. We show that the integrative analysis of genomic and exposomic data has a great potential for understanding genetic and environmental contributions to complex traits and for improving phenotypic prediction accuracy, and thus holds a great promise for future clinical practice.

Methods

Ethics statement

We used data from the UK Biobank (http://www.ukbiobank.ac.uk/) for our analyses. The UK Biobank’s scientific protocol has been reviewed and approved by the North West Multi-centre Research Ethics Committee (MREC), National Information Governance Board for Health & Social Care (NIGB), and Community Health Index Advisory Group (CHIAG). UK Biobank has obtained informed consent from all participants. Our access to the UK Biobank data was under the reference number 14575. The research ethics approval of the current study was obtained from the University of South Australia Human Research Ethics Committee. All methods were performed in accordance with the relevant guidelines and regulations.

Genotype data

The UK Biobank contains health-related data from ~ 500,000 participants aged between 40 and 69, who were recruited throughout the UK between 2006 and 2010[41]. Prior to data analysis, we applied stringent quality control to exclude unreliable genotypic data. We filtered SNPs with an INFO score (used to indicate the quality of genotype imputation) < 0.6, a MAF < 0.01, a Hardy–Weinberg equilibrium p-value < 1e−4, or a call rate < 0.95. We then selected HapMap phase III SNPs, which are known to yield reliable estimates of SNP-based heritability[42-44], for downstream analyses. We filtered individuals who had a genotype-missing rate > 0.05, were non-white British ancestry, or had the first or second ancestry principal components outside six standard deviations of the population mean. We also applied quality control on the degree of relatedness between individuals by excluding one of any pair of individuals with a genomic relationship > 0.025. From the remaining individuals, we selected those who were included in both the first and second release of UK Biobank genotype data. Eventually, 288,837 individuals and 1,133,273 SNPs passed the quality control of genotype data. Among these, 38,921 individuals had no missing data for any of the exposomic variables used in the present study, which were used in the main analysis. Depending on the missingness of the main phenotypic data, sample size varies across traits (see Table 1 for N).

Phenotype data

We chose eleven UK Biobank traits available to us that have a heritability estimate (by an independent open source; https://nealelab.github.io/UKBB_ldsc/) greater than 0.05. These traits are standing height, sitting height, body mass index, heel bone mineral density, fluid intelligence, weight, waist circumference, hip circumference, waist-to-hip ratio, diastolic blood pressure and years of education.

Prior to model fitting, phenotypic data were prepared using R (v3.4.3) in three sequential steps: (1) adjustment for confounders such as age, sex, birth year, social economic status (by Townsend Deprivation Index), population structure (by the first ten principal components of the genomic relationship matrix estimated using PLINK v1.9), assessment centre, and genotype batch using linear regression; (2) standardization; and (3) removal of data points outside  ± 3 standard deviations from the mean.

Exposomic variables

We deliberately selected lifestyle-related variables that are known to affect some of the selected traits to construct the exposome in this study. These variables include smoking, alcohol intake, physical activity, and dietary composition. Details of these variables are listed in Supplementary Table 6. Our aim here is not to cover a comprehensive set of exposomic variables, but to demonstrate the potential use of the proposed integrative analysis of genomic and exposomic data for partitioning phenotypic variance and phenotypic prediction.

Statistical models

We used multiple random-effects LMMs to simultaneously model the effects of the genome and the exposome (model equation ii in Table 2). Genome-exposome correlation was also modelled (model equation iii in Table 2) where the kernel matrix for genome-exposome correlation was explicitly constructed using Cholesky decompositions of g and e[22]. In these models, a genomic relationship matrix (G) was constructed using an n x m1 genotype coefficient matrix (A) as G = , where n is the number of participants and m1 is the number of SNPs. Similarly, an exposomic relationship matrix (E) was estimated using an n x m2 exposomic variable matrix (B) as E= where m2 is the number of exposomic variables (Table 2). These relationship matrices were used to estimate the additive effects of the genome and the exposome. In addition, interaction effects, including gxe, exe and exc, were also considered in these multiple random-effects models (Table 2). The kernel matrices for the interaction terms were derived by the Hadamard product of g and e or e and e (model equations iv, v and vi in Table 2). Reaction norm model[31,32] was used to estimate exc. (model equation v in Table 2). All variance components were estimated using restricted maximum likelihood (REML)[3].

Simulations

Using simulations, we identified two conditions that can cause biased variance estimates of additive effects of exposomic variables and exe interactions, which are correlations between exposomic variables and skewed distributions of exposomic variables. To show the impact of the correlation structure of exposomic variables on variance estimates of exposomic effects, we simulated, for 5,000 individuals, a set of ten orthogonal exposomic variables and another set of ten correlated exposomic variables, each from a multivariate normal distribution. Based on each set of exposomic variables, we then simulated phenotypes using the model y = e + exe + ε with , , and being set to 0.4, 0.5, and 0.1 respectively. The simulated exe effect was based on all possible interaction terms between exposomic variables (as specified in model v of Table 2). The simulation was repeated 100 times, resulting in 100 replicates, each with phenotypes for 5,000 individuals. For each replicate, we fitted the model y = e + exe + ε and averaged variance component estimates across replicates.

Variance estimates of phenotypes that were simulated using correlated and uncorrelated exposomic variables are summarized in Supplementary Table 5. When exposomic variables are orthogonal, all variance-component estimates are unbiased. By contrast, when exposomic variables are correlated, is over estimated, although the estimate of is unbiased. To remedy the effect of correlated exposomic variables on estimate, we used all principal components (PCs) of the correlated exposomic variables to construct the kernel matrix for estimating , and used all pair-wise interaction terms of these PCs to construct the kernel matrix for estimating . Importantly, while retaining all information of the original exposomic variables, the PCs are orthogonal to each other (Jolliffe, 1982). We found that variance estimation based on the PCs of the correlated exposomic variables are unbiased (last column of Supplementary Table 5).

To show the impact of skewness of the distributions of exposomic variables on variance component estimation, we repeated the above simulations using 10 exposomic variables from the UK biobank with skewed distributions. We also noted that these exposomic variables are correlated. Results are presented in Supplementary Table 7. Estimation based on these exposomic variables is biased for both and . Using the PCs of these exposomic variables did not completely eliminate the bias, indicating that skewness of the distributions of exposomic variables affects variance estimation independently from the correlation structure of exposomic variables. As a remedy, we reduced the skewness by removing outliers outside 3 standard deviations from the mean. We found that after this quality control procedure the estimate of became unbiased; but the estimate of remained biased. These results indicate that the estimation of is sensitive to the correlation structure of exposomic variables, while the estimation of is sensitive to the skewness of the distributions of exposomic variables. When using all principal components of the skewness-corrected exposomic variables, all variance estimates became unbiased. Taken together, to avoid biased variance estimation of exposomic effects, it is necessary to (1) conduct quality control on the exposomic variables where values outside 3 standard deviations from the mean should be removed; and (2) transform quality-controlled exposomic variables using a principal component analysis.

We also tested the effect of the correlation structure of exposomic variables on estimate. To do so, we simulated phenotypes based on ten correlated (but quality-controlled) exposomic variables for 5,000 individuals using the model y = g + e + gxe + ε with , , , and being set to 0.3, 0.3, 0.3, and 0.1 respectively. The simulated genetic effect was based on 10 K SNPs that were selected randomly from the 1.1 M Hapmap3 SNPs used for real data analyses, and the simulated gxe effect was based on all possible pairwise interactions between causal SNPs and exposomic variables (as specified in model iv of Table 2). We repeated the simulation 100 times, resulting in 100 replicates. For each replicated, we fitted the model y = g + e + gxe + ε to the genetic data (i.e., 1.1 M Hadmap3 SNPs) and the exposomic data selected for the simulation to estimate variance components, and averaged variance estimates across replicates.

Similar to , the estimation of is affected by the correlation structure of exposomic variables. As shown in Supplementary Table 8, all variance components are biased when the estimation is based on the correlated exposomic variables. Using PCs of the correlated variables corrected the bias for and (see ‘pc1’ in Supplementary Table 8). This observation holds for simulations under a different parameter setting (see Supplementary Table 9) and for simulations based on 10 correlated exposomic variables whose values were simulated from a multivariate normal distribution with a variance–covariance matrix that contains non-zero off-diagonal entries (see Supplementary Table 10 for results).

In Results, we reported for the real data that accounting for significant gxe interactions did not lead to phenotypic prediction accuracy improvements. We hypothesized that the power of phenotypic prediction based on gxe interactions is low; subsequently we used simulations to investigate the power of gxe-based phenotypic prediction. Specifically, we examined, for a sample of 10,000 individuals, of which 80% serves as the training set and 20% as the target set, the extent to which varying effect size of gxe interactions can improve phenotypic prediction accuracy. To do so, we simulated phenotypes using the model y = g + e + gxe + ε with set to 0.2, 0.05, and 0.025, respectively. Each setting has 100 replicates, and each replicate contains phenotypes of 10,000 individuals. We randomly divided each replicate into a training set (n = 8000) and a target set (n = 2000) and subsequently computed the phenotypic prediction accuracy of two estimation models, y = g + e + ε (i.e., null model) and y = g + e + gxe + ε (i.e., full model) for each replicate.

Supplementary Fig. 2 presents the prediction accuracies of the two models by simulation setting (2a) and changes in prediction accuracy from the null model to the full model (2b). Despite the presence of genuine gxe interactions, little prediction accuracy is gained from accounting for these interactions, and this observation holds even under the setting with the largest gxe interactions (i.e., = 0.2). This observation aligns with our results from real data analyses (Supplementary Fig. 1) and indicates that the power of phenotypic predictions based on gxe interactions is low.

Principal component-based transformed variables for E

If the degree of correlation among variables is high, it can cause biased estimates when the variables are fitted in a model, i.e. multicollinearity problem. Such bias is also problematic when using correlated exposomic variables to construct E to be fitted in an LMM to estimate the proportion of the variance explained by the variables (R2 = when phenotypes are standardised with mean zero and variance one). The R2 can also be obtained from a linear model, i.e., the coefficients of determination. For problematically correlated variables, principal component regression has been introduced[45].

A linear model can be written aswhere y is a n vector of phenotypes, W is a column-standardised n x m matrix containing correlated exposomic variables, β is their effects and ε is a vector of residuals.

When exposomic variables in W are highly correlated, estimated exposomic effects (β-hat) are inflated due to multicollinearity problem.

From the singular value decomposition, W can be expressed aswhere U is a matrix whose columns contain the left singular vectors of W, D is a diagonal matrix having a vector containing the singular values of W and V is a unitary matrix (i.e. VVt = I[45]) whose columns contain the right singular vectors of W.

V can be also obtained from the eigen decomposition of the covariance matrix of the variables, i.e. .

The principal component regression approach[45] proposes to transform W to a column-orthogonal matrix, Ω, multiplied by V, which can be written as

Now, we can replace W with Ω in the model as

The estimated regression coefficients from the model (2) are , where V is a unitary matrix such that VVt = I (identity matrix) can be cancelled out[45].

Therefore, R2 values obtained from models (1) and (2) are equivalent as.

However, Eq. (2) can avoid a collinearity issue among the variables. Therefore, model (2) can be extended to a linear mixed model, i.e. the covariance structure can be constructed based on Ω, i.e. ΩΩt/m where Ω is column-standardised.

Suppose a LMM of the formwhere y is a vector of phenotypes for n individuals; W is a n x m2 matrix that contains m exposomic variables; β is a vector of random exposomic effects, each assumed normal with mean zero and variance ; and ε is a vector of residuals, each assumed normal with mean zero and variance .

Under this model, phenotypic variance is partitioned aswhere I is the n × n identify matrix.

When exposomic variables are highly correlated, a transformed W, denoted as Ω, should be used, to avoid biased .

In a similar manner to the linear models (1) and (2), LMM (3) can be rewritten as

Since

Then

Therefore, using column-standardized principal components of exposomic variables as W for Eq. (3) can avoid biased . This is further verified using simulations.

Estimation of exc interactions

We extend the proposed model to a reaction norm model (RNM[31-33,46]) by introducing exc interaction terms, where e is the additive effects of exposomic variables and c is a covariate. Given the significant additive effects found in the above, the interest of fitting RNMs is determine whether these effects vary depending on covariates, which would be evidenced by the presence of significant exc interactions.

While estimates of inform the phenotypic variance explained by the sum of all possible combinations of pairwise interactions between lifestyle-exposomic variables, it may also be of interest to estimate the modulated exposomic effects specific to particular covariates, using the RNM[31-33,46]. The covariates include alcohol intake, smoking, energy intake, physical activity, sex, socio-economic status (indexed by Townsend deprivation index), age and ethnicity measured using the first two ancestry principal components. For each covariate, we fitted the RNM that allows the covariate to interact with exposomic effects and compared the fit of this model with a null model that assumes no exc interactions (see Supplementary Table 3 for p-values). Significant covariates were then included in a subsequent analysis of RNM that fit multiple covariates simultaneously. We reported the total variance of exc interaction effects in Supplementary Table 4.

Five-fold cross-validation

Using fivefold cross validation, we (1) validate significant variance components identified above (Table 1) and (2) evaluate the extent to which the inclusion of these variance components improves phenotypic prediction. For every trait, we randomly split the sample into a discovery set (~ 80%) and a target set (~ 20%) and iterated this process five times in a manner such that target sets did not overlap across iterations (see Fig. 6 for an outline). We derived the prediction accuracy of each model by averaging the Pearson’s correlation coefficients between the observed and predicted phenotypes across target sets; then compared prediction accuracies between models (e.g., y = g + ε vs. y = g + e + ε) to determine phenotypic prediction improvements gained by the inclusion of a given variance component [e.g., var(e)]. For each variance component, we regressed prediction accuracy improvements on estimates of the variance component and declared the variance component valid if the slope differs from zero.

Figure 6

A schematic showing fivefold cross-validation procedures. (i) Randomly assign individuals to 5 groups of an equal size. (ii) Choose one group as the target dataset and the remaining four as the discovery dataset. Iterate the selection process five times in such a way that target datasets do not overlap across iterations. Fit 4 models to each discovery dataset. (iii) For each model, generate the best linear unbiased predictions from discovery datasets and project them onto their corresponding target datasets to derive predicted phenotypes. Compute the phenotypic prediction accuracy for each model by averaging Pearson’s correlation coefficients between the predicted and the observed phenotypes across target datasets.

Theoretical prediction accuracy for quantitative traits

Suppose we predict phenotypes of a quantitative trait (e.g., BMI) with SNP-based heritability h2 using a discovery dataset of N individuals. Following previous theoretical derivations[23,27-30,34], the genomic prediction accuracy based on the model y = g + ε can be written aswhere M1 is the effective number of chromosome segments, which is a function of the effective number of population size and can be estimated using the inverse of the variance of off-diagonal elements of genomic relationships (i.e., G in Table 2) between the discovery and target samples[27-30].

Assuming that phenotypes are standardized to have mean zero and variance one, if the total amount of phenotypic variance explained by the exposome is , Eq. 4 can be adapted to describe the prediction accuracy of the model y = e + ε in the formwhere M2 is analogous to M1 and can be thought of as the effective number of (independent) exposomic variables. Similar to M1, M2 can be estimated using the inverse of the variance of the off-diagonal elements of exposomic relationships (E in Table 2) between the discovery and target samples.

Upon establishing an agreement between expected accuracies, based on Eqs. (4) and (5), and observed accuracies for the 11 traits in this study (Supplementary Fig. 4), we proceeded to the prediction accuracy of the proposed integrative analysis of genomic and exposomic data.

Assuming that the genomic and exposomic effects on phenotypes are uncorrelated, the prediction accuracy of the model y = g + e + ε can be written as

Equation (6) is verified by an agreement between the expected and observed prediction accuracies for the 11 traits in this study (Supplementary Fig. 4).

Theoretical prediction accuracy for disease traits

Considering a disease trait with a population prevalence k, we derived the expected prediction accuracy of the model y = g + e + ε for the disease in terms of the correlation coefficient between the true underlying disease liability and predicted values from the model[23,28,34,47], which can then be converted to an AUC value[23-25].

Similar to rg and re, the expected prediction accuracies for the disease on the liability scale, denoted as (for y = g + ε) and (for y = e + ε), can be computed using previous derivations[23,28,34,47] as the followings.where h2 is the SNP-based heritability on the liability scale, N is the discovery sample size, k is the population prevalence, p is the ratio of cases in the discovery sample, and z is the density at the threshold on the standard normal distribution curve.where is the total amount of variance explained by the exposome on the liability scale (i.e., +). Note when .

As in Eq. (6), we combined and to derive the expected prediction accuracy on the liability scale for the disease, denoted as , under the assumption that the genetic effects and exposomic effects are uncorrelated.

Following a well-established theory[23-25,28] that has been verified by a comprehensive analysis of real data[26], we converted to the area under the receiver operating characteristic curve (AUC) aswhere i (= z/k) is the mean liability for diseased individuals, i2 (= − ik/(1 − k)) is the mean liability for non-diseased individuals, t is the threshold on the normal distribution that truncates the proportion of disease prevalence k and Ф is the cumulative density function of the normal distribution.

To derive the AUC values shown in Fig. 5, we set p = k, M1 to 50,000 and M2 to 28. M1 (50,000) was estimated from the inverse of the variance of genomic relationships (G) between the discovery and target samples[27,29,30]. Similarly, M2 (28) was estimated from the inverse of the variance of exposomic relationships (E) between the discovery and target samples, which agrees with the number of transformed exposomic variables by a principal component analysis in this study (see the correlated exposomic variables section in “Methods”). Note that setting M2 up to 100 would not yield expected prediction accuracies that notably differ from those from setting M2 = 28.

Supplementary Information.