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.


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:// neale lab. 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 Table 1. P-values for estimated variance components of selected traits. Significant estimates are orangecoloured after applying Bonferroni corrected alpha level for each model comparison = 0.05/66 = 7.6E−04. See Table 2 for details of statistical models. www.nature.com/scientificreports/ 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][24][25][26][27][28][29][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.
Exposomic effects on phenotypes. In line with previous estimation (https:// neale lab. 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 lifestyleexposome on phenotypes of all traits (see Fig. 1 for e and Table 1 for p-values under H 0 σ 2 e = 0 ). 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 H 0 σ ge = 0 ; see Supplementary Table 1 for covariance estimates), which was notably different from the genome-transcriptome correlation 22 .
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 H 0 σ gxe|exe = 0 and under H 0 σ exe|gxe = 0).
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 σ 2 e (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 σ 2 exe and σ 2 gxe , 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 G ⊙ E in where y i is the phenotype, μ is the grand mean, a ij is the standardized SNP genotype at locus j, m 1 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 σ 2 g /m 1 , and ε i is the residual assumed to be normal with mean zero and variance σ 2 ε . Where A is a n x m 1 matrix that contains column-standardized genotypes (see the Matrix Notation), and I is the n x n identity matrix (ii) where b k is the kth exposomic variable, m 2 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 σ 2 e /m 2 . To avoid estimation bias due to multicollinearity, b k is transformed using a principal component analysis (see "Methods") Where is a n x m 2 column-orthogonal and columnstandardised matrix that contains the transformed exposomic variables multiplied by their right singular vectors (see 'Principal component-based transformed variables for E' in "Methods") √ 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 σ 2 g×e /Q , and Q is the total number of interaction terms ( Q = m 1 m 2 ) C can be derived using the following pseudo-code with A = a 1 · · · a m1 ; = b 1 · · · b m2 ; C = c 1 · · · c Q , and q = 1, 2 … Q for i = 1 to m 1 { for j = 1 to m 2 { c q = a i ⊙ b j } } σ 2 g G + σ 2 e E + σ 2 g×e Ŵ+σ 2 ε I 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) y i = μ + g i + e i + ∑ x ip θ p P p=1 + ε i e × e i 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 σ 2 e×e /P , and P is the total number of interaction terms (P = m 2 (m 2 + 1)/2). To avoid estimation bias due to multicollinearity, x p is transformed using a principal component analysis (see "Methods") and X can be derived using the following pseudo-code with = b 1 · · · b m2 ; X = x 1 · · · x P , and p = 1, 2 … P for i = 1 to m 2 { for j = i to m 2 { (vi) y i = µ + g i + e i + g × e i + e × e i + ε i y = µ1 n + g + e + g × e + e × e + ε σ 2 g G + σ 2 e E + σ 2 g×e Ŵ + σ 2 e×e �+σ 2 ε I  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.

Continued
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][32][33] . 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 σ 2 e l /m 2 y = µ1 n + g + e + L l=1 e × c l + ε where e × c l is a n × 1 vector that can be derived by e l ⊙ c l , and e l = 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. σ 2 g = phenotypic variance explained by additive effects of the genome; σ 2 e = phenotypic variance explained by additive effects of the exposome; σ 2 exe = phenotypic variance explained by exposome-by-exposome interactions; and σ 2 y = 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 + ε. www.nature.com/scientificreports/ 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).
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][28][29][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][24][25][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  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  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). n t = sample size of training (or discovery) datasets.

Figure 5.
Expected prediction accuracy of the proposed integrative analysis of genetic and exposomic data for disease traits of different prevalence (k) and heritability (h 2 ) at varying levels of total variance explained by the exposome ( σ 2 e.tot ) and sample size of the discovery dataset (N). Diseases are assumed to have a liability of mean zero and variance 1, and both h 2 and σ 2 e.tot 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. www.nature.com/scientificreports/ 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 σ 2 e.tot to denote the total variance in disease liability explained by additive effects of the exposome and exe interactions as a whole.
When setting σ 2 e.tot to 0-that is, using no exposomic information at all-varying the heritability of disease liability h 2 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 h 2 to 0 (i.e., using no genomic information at all), increasing the value of σ 2 e.tot has the same effects on AUC whether using a discovery sample of 50 k or 500 k individuals. Importantly, AUC consistently improves with increasing σ 2 e.tot 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.
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/ hongl ee0707/ 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., σ 2 e that is common to both StructLMM and IGE (Table 3 and Supplmentary Note 1). Importantly, correlated exposomic variables would cause biased estimation of σ 2 e (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 σ 2 e estimates by StructLMM are  Table 5). The other model parameters such as σ 2 exe , σ 2 exc , and σ g,e 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. ukbio bank. 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][43][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). independent open source; https:// neale lab. 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 m 1 genotype coefficient matrix (A) as G = AA t /m 1 , where n is the number of participants and m 1 is the number of SNPs. Similarly, an exposomic relationship matrix (E) was estimated using an n x m 2 exposomic variable matrix (B) as E= t /m 2 where m 2 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 σ 2 e , σ 2 exe , and σ 2 ε 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, σ 2 e is over estimated, although the estimate of σ 2 exe is unbiased. To remedy the effect of correlated exposomic variables on σ 2 e estimate, we used all principal components (PCs) of the correlated exposomic variables to construct the kernel matrix for estimating σ 2 e , and used all pair-wise interaction terms of these PCs to construct the kernel matrix for estimating σ 2 exe . 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 σ 2 e and σ 2 exe . 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 σ 2 exe became unbiased; but the estimate of σ 2 e remained biased. These results indicate that the estimation of σ 2 e is sensitive to the correlation structure of exposomic variables, while the estimation of σ 2 exe 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 σ 2 gxe 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 σ 2 g , σ 2 e , σ 2 gxe , and σ 2 ε 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. www.nature.com/scientificreports/ 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 σ 2 e , the estimation of σ 2 gxe 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 σ 2 e and σ 2 gxe (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 σ 2 gxe 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., σ 2 gxe = 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 (R 2 = σ 2 e when phenotypes are standardised with mean zero and variance one). The R 2 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 as where 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 as where 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. VV t = 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. W t W.
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 where V is a unitary matrix such that VV t = I (identity matrix) can be cancelled out 45 . Therefore, R 2 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.
(1) where y is a vector of phenotypes for n individuals; W is a n x m 2 matrix that contains m exposomic variables; β is a vector of random exposomic effects, each assumed normal with mean zero and variance σ 2 e /m 2 ; and ε is a vector of residuals, each assumed normal with mean zero and variance σ 2 ε . Under this model, phenotypic variance is partitioned as where I is the n × n identify matrix.
When exposomic variables are highly correlated, a transformed W, denoted as Ω, should be used, to avoid biased σ 2 e . In a similar manner to the linear models (1) and (2), LMM (3) can be rewritten as Since V t V = I Then Therefore, using column-standardized principal components of exposomic variables as W for Eq. (3) can avoid biased σ 2 e . This is further verified using simulations.

Estimation of exc interactions.
We extend the proposed model to a reaction norm model (RNM [31][32][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 σ 2 exe 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][32][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.
Theoretical prediction accuracy for quantitative traits. Suppose we predict phenotypes of a quantitative trait (e.g., BMI) with SNP-based heritability h 2 using a discovery dataset of N individuals. Following previous theoretical derivations 23,[27][28][29][30]34 , the genomic prediction accuracy based on the model y = g + ε can be written as where M 1 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][28][29][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 σ 2 e , Eq. 4 can be adapted to describe the prediction accuracy of the model y = e + ε in the form where M 2 is analogous to M 1 and can be thought of as the effective number of (independent) exposomic variables. Similar to M 1 , M 2 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][24][25] .
Similar to r g and r e , the expected prediction accuracies for the disease on the liability scale, denoted as r ′ g (for y = g + ε) and r ′ e (for y = e + ε), can be computed using previous derivations 23,28,34,47 as the followings.
where h 2 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.  (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. www.nature.com/scientificreports/ where σ 2 e.tot is the total amount of variance explained by the exposome on the liability scale (i.e., σ 2 e +σ 2 exe ). Note σ 2 e.tot = σ 2 e when σ 2 exe = 0. As in Eq. (6), we combined r ′ g and r ′ e to derive the expected prediction accuracy on the liability scale for the disease, denoted as r ′ , under the assumption that the genetic effects and exposomic effects are uncorrelated.
Following a well-established theory [23][24][25]28 that has been verified by a comprehensive analysis of real data 26 , we converted r ′ to the area under the receiver operating characteristic curve (AUC) as where i (= z/k) is the mean liability for diseased individuals, i 2 (= − ik/(1 − k)) is the mean liability for nondiseased 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, M 1 to 50,000 and M 2 to 28. M 1 (50,000) was estimated from the inverse of the variance of genomic relationships (G) between the discovery and target samples 27,29,30 . Similarly, M 2 (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 M 2 up to 100 would not yield expected prediction accuracies that notably differ from those from setting M 2 = 28.

Code availability
The source code for MTG2 v2.18 and example code along with related files for fitting IGE model can be accessed without any restrictions from https:// sites. google. com/ site/ hongl ee0707/ mtg2 or from https:// github. com/ hongl ee0707/ IGE.