$$\\widehat{{\\varvec{V}}}=[\\begin{array}{c}\\widehat{{h}_{I}^{2}} \\\\ \\widehat{{\\rho }_{IG}}\\end{array}]\\sim N([\\begin{array}{c}{h}_{I}^{2}\\\\ {\\rho }_{IG}\\end{array}], \\boldsymbol{\\Sigma }),$$<\/p>\n
\n (4)\n <\/p>\n
where \\(\\boldsymbol{\\Sigma }\\) is the true variance\u2013covariance matrix for the two estimates. Here we plug in the empirical variance\u2013covariance matrix \\(\\widehat{\\boldsymbol{\\Sigma }}\\) obtained via delete-block-wise jackknife. In order to indirectly test \\({h}_{I}^{2}\\), we propose a joint modeling test, to test the BV vector \\(\\left[\\begin{array}{c}{h}_{I}^{2}\\\\ {\\rho }_{IG}\\end{array}\\right]=\\left[\\begin{array}{c}0\\\\ 0\\end{array}\\right]\\), which is equivalent to testing the squared Mahalanobis distance estimated by \\(\\widehat{{d}^{2}}={\\widehat{\\boldsymbol{ }{\\varvec{V}}}}^{{\\varvec{T}}}{\\widehat{\\boldsymbol{\\Sigma }}}^{-1}\\widehat{{\\varvec{V}}}\\).<\/p>\n
Fig.\u00a01<\/b> An illustration of the relationship between statistical power and magnitude of \\({h}_{I}^{2}\\) and \\({r}_{IG}\\)<\/p>\n Fig.\u00a02<\/b> Statistical power of different GE testing methods of LDER-GE-based in simulation studies with various parameters<\/p>\n Fig.\u00a03<\/b> UKBB real data analysis results for the 151 E-Y pairs<\/p>\n For simulation studies using real UKBB genotype data, we have a total of 54 different parameter settings, and none of the three methods dominates the other two across all parameter settings (Fig.\u00a02<\/a>). For example, when \\({r}_{IG}\\) is very low or equal to 0, testing only \\({h}_{I}^{2}\\) achieves the highest power among all three methods. On the contrary, as all parameters \\({h}_{g}^{2}\\), \\({h}_{I}^{2}\\), and \\({r}_{IG}\\) increase, testing only \\({\\rho }_{IG}\\) starts to achieve higher power. And joint modeling test achieves a balance of two single-variate test methods and works the best when the true \\({h}_{I}^{2}\\) term is very small but \\({h}_{g}^{2}\\) and \\({r}_{IG}\\) terms are relatively large (Fig.\u00a02<\/a>). We also evaluated the type-I error rate for three tests and concluded that all three methods control the type-I error rate well (Table\u00a01<\/a>), whether they test \\({h}_{I}^{2}\\) directly or indirectly. Since theoretically no optimal test exists for the G\u2009\u00d7\u2009E interaction testing task across all scenarios, we applied our method to UKBB real data analysis (Methods), and the results demonstrated the efficacy of the joint modeling test (Fig.\u00a03<\/a>). Across the 151 E-Y pairs tested, there are more statistically strong (\u2212\u2009log(p value)\u2009>\u200910) signals for \\({\\rho }_{IG}\\) than \\({h}_{I}^{2}\\) (Fig.\u00a03<\/a>A), which can be partially attributed to the larger narrow-sense heritability of the phenotypes. Overall, testing only \\({h}_{I}^{2}\\) identified 35 signals, and testing only \\({\\rho }_{IG}\\) identified 50 signals, while joint modeling test identified 63 signals, covering all signals found by the other two methods except one signal from \\({\\rho }_{IG}\\) (Fig.\u00a03<\/a>B). A full list for all parameter estimates and p values can be found in Additional file 3: Table S1.\n<\/p>\n Table\u00a01 Type-I error rate at 0.05 level for three LDER-GE-based methods to test \\({h}_{I}^{2}\\), each scenario with 500 replications<\/b>Utilize full LD information to estimate G\u2009\u00d7\u2009E genetic covariance more accurately<\/p>\n We have derived the following moment condition for the product of additive genetic effect (GWAS) association vector and G\u2009\u00d7\u2009E interaction effect (GWIS) association vector (Supplementary Note 1):<\/p>\n $$E\\left({{\\varvec{Z}}}_{{\\varvec{G}}}{{\\varvec{Z}}}_{{\\varvec{I}}}^{{\\varvec{T}}}\\right)=\\sqrt{{N}_{G}{N}_{I}}{\\rho }_{IG}{\\varvec{L}}\/M+2{N}_{S}\\left({\\rho }_{IG}+{\\rho }_{01}\\right){\\varvec{R}}\/\\sqrt{{N}_{G}{N}_{I}},$$<\/p>\n \n (5)\n <\/p>\n where \\({{\\varvec{Z}}}_{{\\varvec{G}}}\\) is the Z-score vector for GWAS with sample size \\({N}_{G}\\), and \\({{\\varvec{Z}}}_{{\\varvec{I}}}\\) is the Z-score vector for G\u2009\u00d7\u2009E interaction effect with sample size \\({N}_{I}\\). Here \\({N}_{S}\\) is the sample overlap of two association analysis, \\({\\varvec{R}}\\) is the LD matrix for all variants, and \\({\\varvec{L}}={{\\varvec{R}}}^{{\\varvec{T}}}{\\varvec{R}}\\) is the LD score matrix. We conducted eigen-decomposition on LD matrix as \\({\\varvec{R}}={\\varvec{U}}{\\varvec{D}}{{\\varvec{U}}}^{T}\\), where \\({\\varvec{D}}\\) has diagonal eigenvalues and \\({\\varvec{U}}\\) has eigenvectors and is orthogonal. We transform the original Z-score vectors as \\({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{G}}}^{\\boldsymbol{ }}={{\\varvec{D}}}^{\\left(-1\/2\\right)}{{\\varvec{U}}}^{{\\varvec{T}}}{{\\varvec{Z}}}_{{\\varvec{G}}}\\) and \\({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }}={{\\varvec{D}}}^{\\left(-1\/2\\right)}{{\\varvec{U}}}^{{\\varvec{T}}}{{\\varvec{Z}}}_{{\\varvec{I}}}\\) and obtain the moment condition of the jth diagonal element of the transformed vector product:<\/p>\n $$E{\\left({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{G}}}^{\\boldsymbol{ }}{\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }{\\varvec{T}}}\\right)}_{j,j}=\\sqrt{{N}_{G}{N}_{I}}{\\rho }_{IG}{D}_{j,j}\/M+\\frac{2{N}_{S}\\left({\\rho }_{IG}+{\\rho }_{01}\\right)}{\\sqrt{{N}_{G}{N}_{I}}}.$$<\/p>\n \n (6)\n <\/p>\n The eigen-decomposition and the subsequent Z-vector transformation enable us to exploit more genetic architecture information and result in higher statistical efficiency in estimation.<\/p>\n Formulas (5<\/a>) and (6<\/a>) indicate that sample overlap between GWAS and GWIS does not affect the parameter estimates of \\({{\\varvec{\\rho}}}_{{\\varvec{I}}{\\varvec{G}}}^{ }\\) as \\({{\\varvec{N}}}_{{\\varvec{S}}}\\) appears only in the intercept term. This means that the GWAS and GWIS datasets can be identical, separate, or partially overlapping. In the UKBB application, the same set of individuals was used to compute both \\({{\\varvec{Z}}}_{{\\varvec{G}}}\\) and \\({{\\varvec{Z}}}_{{\\varvec{I}}}\\), consistent with this property.<\/p>\n Through simulations, we show that the empirical standard deviation and root mean square rate of BV-LDER-GE estimate is consistently lower than PIGEON across all simulation scenarios. The average empirical standard deviation across all simulation scenarios of BV-LDER-GE is 23.7% less than PIGEON, approximately equivalent to a sample size increase of 53% (Additional file 3: Table S2). The root mean squared error of BV-LDER-GE is also consistently lower than PIGEON across all simulation scenarios (Additional file 3: Table S2). We selected three simulation scenarios \\({(h}_{I}^{2}{=0.01, h}_{g}^{2}{=0.1,r}_{IG}=0.5),{(h}_{I}^{2}{=0.02, h}_{g}^{2}{=0.2,r}_{IG}=0.5)\\), and \\({(h}_{I}^{2}{=0.03, h}_{g}^{2}{=0.3,r}_{IG}=0.5)\\) to visualize the estimates distribution, statistical efficiency, and unbiasedness (Fig.\u00a04<\/a>). For UKBB real data, BV-LDER-GE detected G\u2009\u00d7\u2009E genetic covariance on 50 E-Y pairs while PIGEON only identified 40. For the 151 E-Y pairs, the average reported standard error of BV-LDER-GE for the 151 E-Y pairs was 30.8% lower than PIGEON (Additional file 3: Table S2), larger but consistent with simulation results.<\/p>\n Fig.\u00a04<\/b> Statistical efficiency comparison of GE genetic covariance estimation for Bivariate-LDER-GE and PIGEON in simulations using real UKBB genotype panel<\/p>\n The authors of PIGEON [7<\/a>] showed that assuming perfectly fitted PRS, the expected value of the coefficient of the PRS-by-E term in a linear regression model equals to the ratio of G\u2009\u00d7\u2009E genetic covariance and the narrow-sense heritability of the phenotype:<\/p>\n $$Y={\\beta }_{1}PRS+{\\beta }_{2}E+{\\beta }_{3}PRS*E+\\upepsilon ,$$<\/p>\n $$E\\left(\\widehat{{\\beta }_{{3}_{\\text{lm}}}}\\right)=\\frac{{\\rho }_{GE}}{{h}_{g}^{2}}.$$<\/p>\n We conducted linear regression analysis on UKBB subjects using fitted PRS from external GWAS summary statistics on the 13 phenotypes and the seven environment covariates to compare the PRS-by-E regression coefficients and G\u2009\u00d7\u2009E genetic covariance estimates using BV-LDER-GE. The PRS fitting was conducted using SDPR [18<\/a>]. The details for fitting PRS, linear regression, and coefficients comparison can be found in Supplementary Note 2. As shown in Fig.\u00a05<\/a>, the estimates are highly correlated (R2\u2009=\u20090.808) and in all E-Y pairs where both methods are significant, the signs of PRS-by-E effects are consistent. This demonstrates the reliability of using additive genetic effect and G\u2009\u00d7\u2009E interaction summary statistics to study the PRS-by-E effect on the phenotype. Note that the estimates of the PRS-by-E regression coefficient is smaller in magnitude than estimated from BV-LDER-GE (Fig.\u00a05<\/a>, Additional file 3: Table S3) due to the attenuation bias resulted from imperfectly fitted PRS [7<\/a>, 19<\/a>] with measurement error. Our results indicate the potential of a more convenient and feasible alternative to evaluate the PRS-by-E effect without fitting PRS [7<\/a>]. On the other hand, the p value of PRS-by-E regression test can be substantially lower than BV-LDER-GE (Additional file 3: Table S3), because BV-LDER-GE only utilizes summary-level statistics instead of individual-level genotype data. Such PRS-by-E effect direction and magnitude analysis may provide biological interpretation on how the environmental and genetic factors shape the outcome together.<\/p>\n Fig.\u00a05<\/b> The PRS-by-E effect coefficient estimates from Bivariate-LDER-GE and linear regression<\/p>\n Some of our G\u2009\u00d7\u2009E interaction effect conclusions are consistent with the literature. For example, type 2 diabetes (T2D) has significant marginal PRS-by-E signals on all seven environmental covariates studies. Among them, the PRS-by-E interaction effect with sex [20<\/a>, 21<\/a>] and smoking [4<\/a>, 22<\/a>] has been reported. Meanwhile, we also report novel interaction discoveries: increased age, being a female, and having higher BMI and drinking frequency magnify the genetic effect of T2D in the European population represented by UKBB. We further ran a conditional PRS-by-E regression including all seven environmental covariates and their PRS-by-E terms in one regression model. Pm2.5 and Townsend deprivation index became nonsignificant while the other five environmental covariates remained significant with multiple testing correction (Additional file 3: Table S4). This demonstrates that multiple environmental exposures and health-related factors interact with genetic predisposition of T2D risk, which might inspire future research direction to the mechanism and intervention regarding T2D and related diseases.<\/p>\n Correct the confounding effects introduced by heritable environmental variable<\/p>\n One important assumption about the model is that the environment variable is either non-heritable or is polygenically heritable in which its additive genetic effects are uncorrelated with the G\u2009\u00d7\u2009E interaction effects on the phenotype. Many other studies consider heritable phenotypes as the interactive environmental variable [7<\/a>, 23<\/a>], such as BMI and alcohol intake [11<\/a>, 12<\/a>], and there is a need to test and correct for the potential confounding effect. PIGEON [7<\/a>] developed a model and framework on this, but utilized only partial LD information through LD score regression [15<\/a>]. And we improve the existing method for more accurate correction using full LD information. We write the heritable E as<\/p>\n $${E}_{i}={\\sum }_{j=1}^{M}{G}_{ji}{\\alpha }_{j}+{\\epsilon }_{2},$$<\/p>\n \n (7)\n <\/p>\n the genetic random effects as<\/p>\n $$[ \\begin{array}{c}{\\beta }_{j}\\\\ {\\gamma }_{j}\\\\ {\\alpha }_{j}\\end{array} ]{\\sim }^{iid} N ( [ \\begin{array}{c}0\\\\ 0\\\\ 0\\end{array} ], [ \\begin{array}{ccc}{h}_{g}^{2}\/M& {\\rho }_{IG}\/M& {\\rho }_{IE}\/M\\\\ {\\rho }_{IG}\/M& {h}_{I}^{2}\/M& {\\rho }_{GE}\/M\\\\ {\\rho }_{IE}\/M& {\\rho }_{GE}\/M& {h}_{E}^{2}\/M\\end{array} ] ),$$<\/p>\n \n (8)\n <\/p>\n and the moment condition for the transformed Z-score product (Supplementary Note 1),<\/p>\n $$E{\\left({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{E}}}^{\\boldsymbol{ }}{\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }{\\varvec{T}}}\\right)}_{j,j}=\\sqrt{{N}_{E}{N}_{I}}\\left({\\rho }_{IE}\\left(1+{h}_{E}^{2}\\right)\\right){D}_{j,j}\/M+{c}_{1},$$<\/p>\n \n (9)\n <\/p>\n $$E{\\left({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }}{\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }{\\varvec{T}}}\\right)}_{j,j}={N}_{I}\\left({h}_{I}^{2}+2{{\\rho }_{IE}}^{2}\\right){D}_{j,j}\/M+{c}_{2},$$<\/p>\n \n (10)\n <\/p>\n $$E{\\left({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{G}}}^{\\boldsymbol{ }}{\\widetilde{{\\varvec{Z}}}}_{{\\varvec{I}}}^{\\boldsymbol{ }{\\varvec{T}}}\\right)}_{j,j}=\\sqrt{{N}_{G}{N}_{I}}{(\\rho }_{IG}+{\\rho }_{GE}{\\rho }_{IE}){D}_{j,j}\/M+{c}_{3},$$<\/p>\n \n (11)\n <\/p>\n where \\({\\widetilde{{\\varvec{Z}}}}_{{\\varvec{E}}}^{\\boldsymbol{ }}={{\\varvec{D}}}^{\\left(-1\/2\\right)}{{\\varvec{U}}}^{{\\varvec{T}}}{{\\varvec{Z}}}_{{\\varvec{E}}}\\) is the transformed Z-score vector for the GWAS on E, with sample size\\({N}_{E}\\), and \\({c}_{s}\\) are the simplified intercept terms. Identifying the slope coefficients of Eqs. (9<\/a>)\u2013 (11<\/a>), we iteratively estimate\\({h}_{E}^{2}, {\\rho }_{IE}, {h}_{I}^{2}, {\\rho }_{GE}\\), and \\({\\rho }_{IG}\\) and use delete-block wise jackknife for inference. We finally use the corrected \\({h}_{I}^{2}\\) and \\({\\rho }_{IG}\\) to conduct the joint modeling test on the genome-level G\u2009\u00d7\u2009E interaction proportion or estimate its magnitude directly. Our simulations showed that both BV-LDER-GE and PIGEON produce unbiased corrected estimation for G\u2009\u00d7\u2009E interaction proportion and G\u2009\u00d7\u2009E genetic covariance, but the estimation efficiency of BV-LDER-GE was 22.2% higher than PIGEON for G\u2009\u00d7\u2009E interaction proportion and 24.7% for G\u2009\u00d7\u2009E genetic covariance (Fig.\u00a06<\/a>, Additional file 3: Table S5).<\/p>\n Fig.\u00a06<\/b> Statistical efficiency comparison of adjusted estimation of GE interaction proportion and GE genetic covariance with heritable environmental covariate confounding effects<\/p>\n We applied two debiased methods to the 151 E-Y pairs of UKBB dataset. Through the 151 tests, BV-LDER-GE identified 30 pairs with G\u2009\u00d7\u2009E interaction proportion bias introduced by the genetic components of the environmental variables, covering all 19 pairs of PIGEON. The confounding test is equivalent to testing term \\({\\rho }_{IE}\\) (Eq.\u00a0(10<\/a>)). After adjusting for potential confounding effects, the BMI-CHO and BMI-LDL pairs failed to pass the joint modeling test of BV-LDER-GE and the rest 61 positive unadjusted BV-LDER-GE joint modeling test remained significant. The details of these 151 adjusted E-Y pairs can be found in Additional file 3: Table S6. As Fig.\u00a07<\/a>A shows, the adjusted \\({h}_{I}^{2}\\) estimates are always smaller than the unadjusted estimates, consistent with the theoretical derivation (Eq.\u00a0(10<\/a>)). And the p values will increase if the confounding effect is non-trivial (Fig.\u00a07<\/a>B). For common practice, we recommend testing the existence of term \\({\\rho }_{IE}\\) first as some social economic status covariates may still contain genetic components [24<\/a>]. If the test is significant then an adjustment of genetic confounding effect should be considered, even if a significant confounding effect may still lead to small bias (Fig.\u00a07<\/a>A).<\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig1_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig2_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig3_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig4_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig5_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig6_HTML.png\")
<\/a><\/p>\n
![\"figure](\"https:\/\/www.europesays.com\/uk\/wp-content\/uploads\/2025\/10\/13059_2025_3815_Fig7_HTML.png\")
<\/a><\/p>\n