MainGenome-wide association studies (GWAS) have successfully identified many genetic variants associated with diseases and disease-related phenotypes, reshaping our understanding of the genetic architecture underlying complex traits. Most GWAS, however, have focused on genetic main effects, implicitly assuming a constant effect size across environmental contexts and overlooking the fact that genetic effects can be context-dependent and modulated by environmental factors1,2. Indeed, both genetic and environmental factors are major contributors to phenotypic variation, and their interactions play a pivotal role in shaping complex human traits3,4. For example, the genetic predisposition to obesity may manifest differently depending on lifestyle factors such as diet and physical activity5. The recent availability of large-scale biobank data has provided an unprecedented opportunity to systematically characterize genotype–environment interactions (GxE) at a scale and resolution previously unattainable. Leveraging these datasets, GxE studies have started to dissect the dynamic interplay between genes and the environment6, uncover context-specific genetic effects7, generate deeper insights into the architecture of complex traits8, build more accurate models for genetic risk prediction by incorporating GxE effects9, and ultimately pave the way for personalized medicine strategies that integrate both genetic and environmental factors10,11.Despite the promise of GxE studies, conducting effective genome-wide GxE association analysis in biobank scale studies presents several important statistical and computational challenges. Two primary strategies have been developed for detecting GxE associations. The first strategy focuses on single-environment GxE tests, using tools such as the linear regression model framework implemented in PLINK12 and the linear mixed model framework implemented in fastGWA-GE13. However, these approaches often lack sufficient statistical power, even in large-scale biobank data, and have sometimes failed to detect well-established GxE loci, such as the FTO locus for body mass index (BMI)13. An alternative strategy involves joint GxE testing across multiple environmental factors, as implemented in tools such as StructLMM14 and GEM15. These multi-environment GxE approaches offer much improved statistical power compared with single-environment GxE tests, as genetic variants often interact with multiple environmental factors to influence phenotypes6,14. However, these approaches have three key limitations. First, their statistical power can be substantially reduced when only a small subset of environments contributes to GxE interactions. Second, they are computationally intensive, restricted to analyzing candidate genomic regions, and lack the scalability required for genome-wide GxE analysis in biobank-scale datasets. Third, these approaches do not account for polygenic effects or polygenic interaction effects arising from genome-wide variants, nor do they model heterogeneous noise variance across environments. As we demonstrate here, both types of polygenic effects explain a substantial proportion of phenotypic variance, and failing to account for them can lead to excessive false positives in GxE association analyses13,16. In addition, ignoring heterogeneous noise can bias estimates of the underlying genetic architecture17.Here we present fast genotype–environment interaction analysis (fastGxE), a scalable and effective genome-wide multi-environment GxE method designed to identify genetic variants that interact with one or multiple environmental factors. fastGxE controls for polygenic effects, polygenic interaction effects and noise heterogeneity, remains robust to the number of interacting environments, and is 32.98–126.49 times faster than existing approaches. We also developed mmSuSiE, an extension of SuSiE18 specifically designed for mixed-model analysis, to identify the environmental factors driving the GxE interactions detected by fastGxE. We illustrate the benefits of fastGxE through extensive simulations and an in-depth analysis of 32 physical traits and 67 blood biomarkers from the UK Biobank.ResultsfastGxE overview and runtime comparisonfastGxE is described in detail in the Methods, with technical specifics provided in the Supplementary Note. In brief, fastGxE is a scalable and effective method designed for genome-wide GxE association analysis (Fig. 1). fastGxE examines one SNP at a time, testing its interactions with multiple environmental factors while controlling for their main effects, polygenic effects, polygenic interaction effects and noise heterogeneity across multiple environmental factors. This enables fastGxE to produce well-calibrated P values while yielding high statistical power. For SNPs identified as candidate GxE signals, we further developed a mixed-model SuSiE method, mmSuSiE, which is an extension of SuSiE18, to identify the specific environmental factors driving the observed GxE interactions.Fig. 1: Schematic overview of fastGxE.Full size imagefastGxE is a scalable and effective method designed for genome-wide GxE association analysis. a, fastGxE handles multiple environmental factors and examines one SNP at a time, decomposing the phenotype into the SNP main effect, environmental main effects and GxE interaction effects, while controlling for polygenic effects, polygenic interaction effects and noise heterogeneity. b, fastGxE evaluates various GxE effect size configurations and combines the resulting P values into a single P value to test whether the SNP interacts with at least one environmental factor. c, By explicitly modeling polygenic background and heteroscedastic noise, fastGxE generates calibrated P values for identifying candidate GxE loci. d, In addition, fastGxE utilizes mmSuSiE, an extension of the SuSiE algorithm, to identify the environmental factors driving the detected GxE interactions and uses the stratified Wald test to visualize and support these interactions. In a and b, \({\boldsymbol{y}}\) is an n-vector of phenotypic measurements; \({{\boldsymbol{x}}}_{{\boldsymbol{G}}}\) is an \(n\)-vector of genotypes for the SNP of focus; \({{\boldsymbol{X}}}_{E}\) is an n by q matrix for \(q\) measured environmental factors; \({\beta }_{G}\) is the genetic main effect size; \({{\boldsymbol{\beta }}}_{E}\) is a q-vector of environment main effect sizes; \({{\boldsymbol{\beta }}}_{{GxE}}\) is a q-vector of interaction effects between the SNP of focus and \(q\) environmental factors; \(\odot\) denotes the Hadamard product; \({\boldsymbol{g}}\) is an n-vector of polygenic effects; \({\boldsymbol{v}}\) is an n-vector of polygenic GxE effects; \({\boldsymbol{u}}\) is an n-vector of multiple-environment NxE effects; \({\boldsymbol{e}}\) is an n-vector of residual errors. In d, the AIE score summarizes the joint contribution of multiple environments to the GxE interaction at a given SNP and represents predicted individual phenotypic values based on GxE interactions across multiple environments, with high scores indicating greater phenotypic values driven by GxE. In c, the dashed diagonal indicates the null expectation, the gray shaded region illustrates the 95% interval expected under the null, and the dashed horizontal line the significance threshold; in d, orange points denote environmental factors with high PIP, and the vertical bars in the stratified-effect plot illustrate error bars around effect estimates. All panels are schematic: the plots in c and d (QQ plot, locus-zoom plot, PIP plot and stratified-effect plot) are illustrative drawings rather than observed or simulated data, and no sample size, replicates, statistical tests or data-derived error bars or intervals are associated with this figure; the shaded band in c and the bars in d are drawn for illustration only. Figure created in BioRender; Ning, C. https://biorender.com/yy0spyc (2026).We first assessed fastGxE’s computational efficiency by comparing it with the single-environment GxE method fastGWA-GE and the multi-environment GxE method StructLMM across varying numbers of CPU cores, number of environmental factors and sample sizes (Fig. 2a–c). As the number of CPU cores increased, the runtime for each method decreased, but fastGxE exhibited a much greater reduction, demonstrating its superior multithreading performance. In addition, as the number of environmental factors increased, fastGxE’s runtime grew at a much slower rate compared with the other methods, supporting its better scalability. fastGxE also completed the analyses faster even when applied to a single environmental factor, despite being specifically designed as a multi-environment GxE method. This computational efficiency makes fastGxE particularly well suited for large biobank datasets. For example, in an analysis of 300,000 individuals with 7.8 million SNPs across 64 environmental factors, fastGxE completed the analysis in just 6.3 h using 16 cores—32.98 times faster than fastGWA-GE and 126.49 times faster than StructLMM, making it feasible for conducting genome-wide GxE analysis. Moreover, fastGxE implements a block-wise parallelization strategy that efficiently utilizes multiple threads across multiple nodes, which further improves computation time by 1.86 and 3.22 times when using 16 and 32 cores, respectively (Supplementary Fig. 1).Fig. 2: Runtime and simulation comparisons of fastGxE, StructLMM and fastGWA-GE.Full size imagea–c, Computing time/runtime across varying numbers of CPU cores, number of environmental factors and sample sizes for a total of 7,833,003 SNPs. For StructLMM, the computation was performed on 5,000 SNPs, with runtime linearly extrapolated to the full SNP set. The baseline setting consisted of 300,000 individuals, 64 environmental factors and 16 CPU cores. We systematically varied the number of CPU cores (a), the number of environmental factors (b) and the sample size (c) to evaluate the computational scalability of each method. All analyses were conducted on a computing platform equipped with dual Intel Xeon Gold 6138 CPUs, each operating at 2.00 GHz and comprising 40 cores. d–f, QQ plots of −log10P values from SNPs on chromosome 22 without GxE effects, under NxE effect levels of 0% (d), 15% (e) and 25% (f), simulated using SNPs from chromosomes 1–21. g–i, Power comparisons at varying levels of GxE interaction effect size (in GxE heritability) (g), different numbers of active interacting environmental factors (h) or varying numbers of masked active environmental factors (i)—those contributing to GxE interactions but omitted from the analysis. In d–i, each simulation replicate included 100,000 individuals, and all P values were based on two-sided tests. For fastGxE, each P value is the Cauchy combination of the score tests under the single-environment GxE models and the variance-component test under the multi-environment model (Methods); for StructLMM, the P value is from its variance-component test of GxE; for fastGWA-GE, the P value is the minimum over environmental factors of the single-environment tests, corrected for the number of environmental factors tested (Methods). In d–f, each setting comprised 100 independent simulation replicates; P values for all chromosome 22 SNPs were pooled across replicates. The dashed diagonal lines indicate the null expectation, and the gray shaded regions indicate the 95% confidence interval under the null. In g–i, power is the proportion of the 100 replicates in which the causal SNP was detected at a 5% family-wise error rate, obtained by Bonferroni correction over one million independent SNPs.Source dataSimulationsWe performed comprehensive simulations to evaluate the performance of fastGxE and compared it with StructLMM and fastGWA-GE. First, we examined null settings to evaluate their type I error control across varying levels of heterogeneous environmental noise (NxE), where \({\rho }_{{NxE}}\) denotes the proportion of phenotypic variance explained by environment-dependent noise heterogeneity across multiple environmental factors (\({\rho }_{{NxE}}=0 \% ,15 \% ,25 \%\)). Quantile–quantile (QQ) plots and genomic inflation factors (λ) at different quantiles showed that fastGxE was well calibrated (Fig. 2d–f and Supplementary Table 1). By contrast, fastGWA-GE produced calibrated P values at smaller quantiles but was markedly conservative at higher quantiles, whereas StructLMM yielded reasonably calibrated P values at low \({\rho }_{{NxE}}\) but became inflated as \({\rho }_{{NxE}}\) increased. Next, we assessed type I error control under varying SNP heritability (\({h}_{{SNP}}^{2}=5 \% ,30 \% ,50 \%\)) and GxE heritability, the latter defined as the proportion of phenotypic variance explained by polygenic interaction effects \(({h}_{{GxE}}^{2}=1 \% ,5 \% ,10 \%\)). Across quantiles, fastGxE maintained robust type I error control. For StructLMM, increasing \({h}_{{SNP}}^{2}\) had little effect on type I error, whereas increasing \({h}_{{GxE}}^{2}\) led to inflation (Supplementary Fig. 2 and Supplementary Table 1). In addition, we considered a proxy scenario in which the true environment and the causal SNP exerted only main effects, without interaction, but analyses were conducted using a proxy environment and a proxy SNP. Our results show that fastGxE did not produce false-positive signals in the vast majority of settings (Supplementary Fig. 3). Noticeable inflation occurred only under a specific combination of conditions: low minor allele frequency (MAF; 0.01), high SNP heritability (≥1%) and moderate linkage disequilibrium (LD; r2 = 0.3, 0.5, 0.7) between the causal and proxy SNPs. This combination of conditions appears uncommon in practice (see the ‘GxE association analysis for 32 physical traits’ and ‘GxE association analysis for 67 blood biomarkers’ sections below).Next, we compared the power of the methods at a 5% family-wise error rate, focusing on the influence of (i) GxE interaction effect sizes, (ii) the number of interacting environmental factors and (iii) the number of masked environmental factors—that is, factors that contribute to GxE interactions but are omitted from the analysis. Across all scenarios, fastGxE consistently achieved higher power than fastGWA-GE and StructLMM (Fig. 2g–i). As GxE interaction effect sizes decreased, power declined for all methods; however, fastGxE maintained a clear advantage (Fig. 2g). As the number of interacting environmental factors increased, power for fastGWA-GE declined rapidly, whereas the multi-environment methods were more stable, with fastGxE attaining the highest power (Fig. 2h). Under the masked-environment settings, power decreased for all methods as more environmental factors were masked; yet fastGxE remained consistently more powerful than the other two (Fig. 2i). In addition, fastGxE was robust to the inclusion of irrelevant environmental factors in the analysis (Supplementary Fig. 4).Additional simulations incorporating 96 SNP annotations from BaselineLD v2.2 (ref. 19) showed that fastGxE remained well calibrated and retained higher power under more complex genetic architectures (Supplementary Fig. 5). We also assessed whether prefiltering SNPs based on main genetic effects improves power for GxE analysis and found that this strategy was beneficial only when the genetic main effects were sufficiently large or strongly correlated with the GxE effects; otherwise, prefiltering reduced power (Supplementary Fig. 6). We further assessed sample size requirements using simulations guided by real data (Supplementary Note) and found that, to achieve >80% power, ~100,000 individuals are needed to detect large effects (explaining >0.12% of phenotypic variance), ~300,000 for moderate effects (0.04–0.12%), ~800,000 for small effects (0.02–0.04%) and >1,000,000 for very small effects (