A scalable and robust variance components method reveals insights into the architecture of gene-environment interactions underlying complex traits

Abstract

Understanding the contribution of gene-environment interactions (GxE) to complex trait variation can provide insights into disease mechanisms, explain sources of heritability, and improve genetic risk prediction. While large biobanks with genetic and deep phenotypic data hold promise for obtaining novel insights into GxE, our understanding of GxE architecture in complex traits remains limited. We introduce a method to estimate the proportion of trait variance explained by GxE (GxE heritability) and additive genetic effects (additive heritability) across the genome and within specific genomic annotations. We show that our method is accurate in simulations and computationally efficient for biobank-scale datasets. We applied our method to common array SNPs (MAF ≥1%), fifty quantitative traits, and four environmental variables (smoking, sex, age, and statin usage) in unrelated white British individuals in the UK Biobank. We found 68 trait-E pairs with significant genome-wide GxE heritability (p<0.05/200) with a ratio of GxE to additive heritability of ≈6.8% on average. Analyzing ≈8 million imputed SNPs (MAF ≥0.1%), we documented an approximate 28% increase in genome-wide GxE heritability compared to array SNPs. We partitioned GxE heritability across minor allele frequency (MAF) and local linkage disequilibrium (LD) values, revealing that, like additive allelic effects, GxE allelic effects tend to increase with decreasing MAF and LD. Analyzing GxE heritability near genes highly expressed in specific tissues, we find significant brain-specific enrichment for body mass index (BMI) and basal metabolic rate in the context of smoking and adipose-specific enrichment for waist-hip ratio (WHR) in the context of sex.

Keywords: UK Biobank; complex traits; gene-context interaction; gene-drug interaction; gene-environment interaction; genetic architecture of gene-environment interactions; noise heterogeneity; patitioning GxE heritability; scalable variance component analysis.

Figure 1

Calibration and power of GENIE in simulations ($N=291,273$ unrelated individuals, $M=454,207$ SNPs) (A) Q-Q plot of p values (of a test of the null hypothesis of zero GxE heritability) when GENIE is applied to phenotypes simulated in the absence of GxE effects. Each panel contains 100 replicates of phenotypes simulated with additive heritability $h_g^2=0.25$ and varying proportions of causal variants. The causal ratios are the same for the G and GxE components ($10\%$), and the causal SNPs for the GxE component are independently sampled to those for the additive genetic component. Across all architectures, the mean of $P($rejection at $p<t)$ is $7.5\%$ and $0\%$ for $t=0.05$ and $t=\frac{0.05}{200}$, respectively ($7.5\%$ is not significantly different from the nominal rate of $5\%$). (B) The power of GENIE across genetic architectures as a function of GxE heritability. We report power for p value thresholds of $t\in \left\{0.05,\frac{0.05}{200}\right\}$. (C) The accuracy of $h_{gxe}^2$ estimates obtained by GENIE. Across all simulations, statin usage in UKB was used as the environmental variable.

Figure 2

Comparisons of false positive rates with existing methods with the presence of noise heterogeneity False positive rates of tests for GxE heritability across GENIE, MEMMA, and MonsterLM using (A) continuous and (B) discrete environment exposures. We performed simulations with no GxE heritability but with varying magnitudes of the variance of the NxE effect. We computed the false positive rate as the fraction of rejections (p value of a test of the null hypothesis of zero GxE heritability $<0.05$) over 100 replicates of phenotypes. The phenotypes were simulated from $N=40,000$ individuals and $M=223,591$ SNPs filtered from $M=454,207$ SNPs with the genotype QC steps in MonsterLM: SNPs that failed the Hardy-Weinberg test at the significance threshold ${10}^{-10}$ were excluded, and highly correlated SNPs with LD $r^2>0.9$ and SNPs with MAF $<0.05$ were removed. Error bars correspond to the estimated 95% CI of the rejection rate.

Figure 3

Estimation of G and GxE heritability in six simulated scenarios We investigated the performance of GENIE in estimating G and GxE heritability under six simulated scenarios. (1) Correlated Y: the phenotypes were correlated with the continuous environment exposure, with Pearson’s correlation $r=0.5$; (2) heritable E: the environment exposure E was simulated from the same set of genotype data as in the phenotype simulation, with an additive genetic heritability of 0.1; (3) same causal SNPs: additive genetic causal SNPs completely overlap with GxE causal SNPs; (4) same causal SNPs for additive and heritable E: additive genetic causal SNPs completely overlap with the causal SNPs explaining heritability in E, where E is the same as in scenario (2); (5) collider bias: the phenotype Y and environment exposure E are correlated through an unobserved confounder; we simulated a heritable environment variable with a genetic heritability of 0.1. The phenotypes were then generated to have a Pearson’s correlation $r=0.2$ with the heritable E. We assumed that the correlation was due to an unobserved confounder. (6) Heavy-tailed noise: we drew the environment noise component from the Student’s t-distribution with degrees of freedom = 4. In all scenarios, we simulated 100 replicates of phenotypes with NxE and varying magnitude of GxE effects across $N=291,273$ individuals genotyped at $454,207$ SNPs. The ground truth GxE heritability was 0, 0.04, and 0.1, with corresponding NxE variance of 0.04, 0.04, and 0.1. The additive genetic heritability was fixed at 0.25. The x and y axes denote the true GxE heritability and the estimated G and GxE heritability. Points and error bars represent the mean and estimated $95\%$ CI, respectively. Across all simulations where there is no GxE, the mean of P(rejection at $p<t$) are $5.5\%$ and $0\%$ for $t=0.05$ and $t=0.05/200$, respectively ($5.5\%$ is not significantly different from the nominal rate of $5\%$).

Figure 4

Effect of noise heterogeneity (NxE) on estimates of heritability associated with GxSmoking across 50 quantitative phenotypes in UKB Model G + GxE refers to a model with additive and gene-by-environment interaction components where the environmental variable is smoking status. Model G + GxE + NxE refers to a model with additive, gene-by-environment interaction, and noise heterogeneity (noise-by-environment interaction) components. (A) We ran GENIE under G + GxE and G + GxE + NxE models to assess the effect of fitting an NxE component on the additive and GxE heritability estimates. (B) Comparison of GxE heritability estimates obtained from GENIE under a G + GxE + NxE model (x axis) to a G + GxE model (y axis). Black error bars mark $\pm$ standard errors centered on the estimated GxE heritability. The color of the dots indicates whether estimates of GxE heritability are significant under each model. (C) We performed permutation analyses by randomly shuffling the genotypes while preserving the trait-E relationship and applied GENIE in each setting under G + GxE and G + GxE + NxE models. We report the fraction of rejections P(p value of a test of the null hypothesis of zero GxE heritability $<\frac{0.05}{200}$ that accounts for the number of phenotypes tested) over 50 UKB phenotypes.

Figure 5

Estimates of GxE heritability across phenotypes in UKB Estimates of (A) GxSmoking, (B) GxSex, (C) GxAge, and (D) GxStatin heritability across 50 UKB phenotypes. We applied GENIE to $N=291,273$ unrelated white British individuals and $M=454,207$ array SNPs (MAF $\ge 1\%$). Our model includes the environmental variable as a fixed effect and accounts for noise heterogeneity. The environmental variable is standardized in these analyses. Error bars mark $\pm 2$ standard errors centered on the point estimates. The asterisk and double asterisk correspond to the nominal $p<0.05$ and $p<0.05/200$, respectively.

Figure 6

Estimates of the ratio of GxE to additive heritability across phenotypes in UKB Estimates of the ratio of (A) GxSmoking, (B) GxSex, (C) GxAge, and (D) GxStatin to additive heritability across 50 UKB phenotypes. Error bars mark $\pm 2$ standard errors centered on the point estimates. The asterisk and double asterisk correspond to the nominal $p<0.05$ and $p<0.05/200$, respectively.

Figure 7

Per-allele squared GxE and additive effect sizes as a function of MAF and LD (A) The squared per-allele GxE effect size for four selected pairs of trait and environment (trait-E pairs). (B) The squared per-allele additive effect size for the same trait-E pairs. The x axis corresponds to MAF-LD annotations where annotation $i.j$ includes SNPs in MAF bin i and LD quartile j where MAF bin 1 and MAF bin 2 correspond to SNPs with MAF $\le 5\%$ and MAF $>5\%$, respectively, while the first quartile of LD scores correspond to SNPs with the lowest LD scores respectively). The y axis shows the per-allele GxE (or additive) effect size squared defined as $\frac{h_k^2}{2M_k{f}_k\left(1-f_k\right)}$ where $h_k^2$ is the GxE (or additive) heritability attributed to bin k, $M_k$ is the number of SNPs in bin k, and $f_k$ is the mean MAF in bin k. Error bars mark $\pm 2$ standard errors centered on the estimated effect sizes.

Figure 8

Partitioning GxE heritability across 53 tissue-specific genes We plot $-lo{g}_{10}\left(p\right)$ where p is the corresponding p value of the tissue-specific GxE enrichment defined as $\frac{h_{gxe,tissue}^2/h_{gxe,total}^2}{M_{tissue}/M_{total}}$. For every tissue-specific annotation, we use GENIE to test whether this annotation is significantly enriched for per-SNP heritability, conditional on 28 functional annotations that are part of the baseline LDSC annotations. The dashed and solid lines correspond to the nominal $p<0.05$ and FDR $<0.1$ threshold, respectively. We have labeled two tissues with the most significant p values for each figure.