SharePro for GxE analysisSharePro (Fig. 1) takes exposure-stratified GWAS summary statistics as input. By integrating statistical evidence from subpopulations with different exposure statuses, SharePro can effectively identify causal variants that are either shared or specific to different exposure categories, reducing the multiple testing burden in GxE analysis.Fig. 1: SharePro for GxE analysis.We showcase three settings of effect heterogeneity in a population with three environmental exposure statuses. In setting 1, with similar effect sizes across three exposure categories, this genetic effect can be well-characterized through combined GWAS and fine-mapping of the combined GWAS summary statistics. The stratified approaches might not be well-powered due the relatively small sample sizes in each category. Nevertheless, with SharePro, we can fine-map this genetic effect from exposure stratified GWAS summary statistics and assess its effect heterogeneity. In setting 2, with moderate effect size discrepancy, this genetic effect may be characterized through the combined approaches but the joint approaches will have a higher power. Classical GxE analysis can be underpowered due to a high multiple testing burden. With SharePro, we can first fine-map this genetic effect from exposure stratified GWAS summary statistics and then characterize its effect heterogeneity. With reduced multiple testing burden, we can effectively detect GxE in this setting. In setting 3, the combined approaches can no longer adequately characterize the genetic effect as two exposure categories demonstrates opposite effect sizes. With SharePro, we can accurately fine-map this genetic effect and detect its effect heterogeneity. GWAS, genome-wide association study; GxE, gene-environment interaction.Specifically, similar to the original SharePro method for colocalization analysis12, we use effect groups to align causal signal across exposure stratified GWAS summary statistics, thereby accounting for the uncertainty associated with identifying causal variants due to LD. We assume that there are up to K causal signals for a phenotype y across subpopulations with different exposure categories within a locus containing G variants. For the kth causal signal, we use an effect group indicator sk to represent the causal variant and other variants that are highly correlated with it:$${{{{\bf{s}}}}}_{k} \sim {{{\bf{Multinomial}}}}\left(1,{{{{\bf{1}}}}}_{G\times 1}\times \frac{1}{G}\right)$$We use cke to represent the causal status of the kth effect group in the subpopulation with the eth exposure category:$${c}_{ke} \sim {{{\bf{Bernoulli}}}}(\sigma )$$Additionally, we use βke to represent the effect size of the kth effect group in the eth exposure category:$${\beta }_{ke} \sim {{{\mathcal{N}}}}(0,{\tau }_{\beta }^{-1})$$With the genotype matrix and the phenotype in the subpopulation with the eth exposure category denoted as Xe and ye, we have:$${{{{\bf{y}}}}}_{e} \sim {{{\mathcal{N}}}}\left({{{{\bf{X}}}}}_{e}{\sum}_{k}{{{{\bf{s}}}}}_{k}{\beta }_{ke}{c}_{ke},{\tau }_{y}^{-1}{{{\bf{I}}}}\right)$$Here, τβ, τy and σ are hyperparameters in the model. We use similar strategies for estimating hyperparameters as in our previous work11,12 detailed in Supplementary Notes. Using an efficient variational inference algorithm adapted for GWAS summary statistics11,12,13,14, we can obtain variant-level fine-mapping results (posterior inclusion probabilities, PIP) and identify effect groups for GxE analysis as detailed in Supplementary Notes.Accounting for effect heterogeneity improved power for fine-mappingTo examine the benefits of accounting for effect heterogeneity in fine-mapping, we conducted simulation studies with varying sample sizes (Ne and Nu), effect sizes (βe and βu) and number of causal variants (KC). As expected, in the presence of simulated effect heterogeneity (βe ≠ βu), the 2-df association test, which accounts for effect heterogeneity, exhibited a higher AUPRC compared to the 1-df association test. SharePro fine-mapping accounted for both effect heterogeneity and LD, thereby achieving the highest AUPRC. For example, with three causal variants having an opposite effect direction and equal effect size in the exposed and the unexposed category of equal sample size (Ne = Nu  = 25, 000; βe = 0.05; βu  = −0.05 and KC = 3), the 1-df association test and the combined fine-mapping of the main effect with SparsePro and SuSiE both had an AUPRC  < 0.01, close to the AUPRC of a random predictor (Fig. 2 and Supplementary Data 1), because variant-phenotype associations in the overall population were not detectable. In contrast, the 2-df association test achieved a median AUPRC of 0.34, while SharePro fine-mapping achieved a median AUPRC of 0.95 (Fig. 2 and Supplementary Data 1).Fig. 2: Accounting for effect heterogeneity improved power for fine-mapping.Area under precision-recall curve (AUPRC) in simulation settings with 3 causal variants per locus (KC = 3) is displayed. Rows represent different effect size discrepancies with βe and βu corresponding to effect size in the exposed and the unexposed category, respectively. Columns represent different sample size settings with Ne and Nu corresponding to sample size in the exposed and the unexposed category, respectively. 1-degree of freedom (1-df) and 2-df tests were implemented in GEM. SuSiE and SparsePro were applied to the combined GWAS of the main effect (i.e. 1-df tests). Each dot indicates the AUPRC based on one of the three randomly selected loci, aggregating 50 simulation replicates under the corresponding simulation setting. Detailed statistics are provided in Supplementary Data 1.As the opposing effect (βu) weakened and the combined effect in the overall population \(\left(\frac{{N}_{e}{\beta }_{e}+{N}_{u}{\beta }_{u}}{{N}_{e}+{N}_{u}}\right)\) increased, the performance differences between methods using the joint approaches and those employing the combined approaches attenuated. For example, with causal variants having a effect size discrepancy of 0.07 (βe = 0.05; βu = −0.02; Ne = Nu = 25, 000 and KC = 3), SharePro fine-mapping had a median AUPRC of 0.92 while combined fine-mapping with SuSiE and SparsePro had a median AUPRC of 0.22 and 0.27, respectively (Fig. 2 and Supplementary Data 1). Meanwhile, 2-df association test achieved a median AUPRC of 0.32 and 1-df association test only had a median AUPRC of 0.09 (Fig. 2 and Supplementary Data 1). When the effect size difference decreased to 0.05 (βe = 0.05; βu = 0.00), SharePro fine-mapping maintained a median AUPRC of 0.91 and SuSiE and SparsePro achieved a median AUPRC of 0.77 and 0.80 (Fig. 2 and Supplementary Data 1). Furthermore, as the effect heterogeneity decreased with reduced sample size of the unexposed category, the performance differences also attenuated. For example, with three causal variants having an opposite effect in the exposed category and the unexposed category of unequal sample sizes (Ne  = 25,000; Nu = 10,000; βe  =  0.05; βu  = − 0.05 and KC  =  3), the median AUPRC of SharePro fine-mapping was 0.93 while the median AUPRC of SuSiE and SparsePro were 0.40 and 0.39 (Fig. 2 and Supplementary Data 1). With a smaller sample size for the unexposed category (Nu = 5000), the median AUPRC of SharePro fine-mapping was 0.92 while the median AUPRC of SuSiE and SparsePro were 0.81 and 0.82 (Fig. 2 and Supplementary Data 1).These trends were consistent when there was one or two causal variants (Supplementary Figs. 1–2 and Supplementary Data 1). Meanwhile, SharePro achieved the highest AUROC close to 1 across all simulation settings while the other methods could have a low AUROC in the presence of strong effect heterogeneity (e.g. Ne  = Nu = 25000; βe  = 0.05; βu  = − 0.05; Supplementary Figs 3–5 and Supplementary Data 2).SharePro demonstrated improved power in GxE analysisWe next examined the utility of SharePro in GxE analysis and compared it to LD clumping, COJO and variant-level GxE analysis. First, without effect heterogeneity (βe = βu = 0.05), all methods provided calibrated p-values with well-controlled type I error in different simulation settings and under various significance thresholds (Supplementary Data 3). With addtional 50,000 replications, the SharePro effect-group level GxE p-value maintained calibration at lower significance levels (Supplementary Fig. 6).Importantly, SharePro consistently achieved the highest power with the lowest empirical FDR when there is effect heterogeneity. For example, when the effect size discrepancy between the exposed and the unexposed (βe  −  βu) was large (βe = 0.05, βu = − 0.05 and KC  =  3), at a p-value cutoff of 0.01, all methods achieved a similar power close to 1 (Fig. 3 and Supplementary Data 4). However, multiple testing corrected variant-level GEM p-value still had an empirical FDR of 0.98 (Fig. 3 and Supplementary Data 4). In contrast, the SharePro effect-group level GxE p-value had the lowest empirical FDR of 0.01 while the SharePro GxE p-value mapped to top variants in effect groups (SharePro variant) had an empirical FDR of 0.13 (Fig. 3 and Supplementary Data 4).Fig. 3: SharePro demonstrates improved power in GxE analysis.Power and empirical false discovery rate (FDR) in GxE analysis are derived with a Benjamini-Hochberg-corrected two-sided p-value cutoff of 0.01 with 3 causal variants per locus (KC = 3). Rows represent different effect size discrepancies with βe and βu corresponding to effect size in the exposed and the unexposed category, respectively. Columns represent different sample size settings with Ne and Nu corresponding to sample size in the exposed and the unexposed category, respectively. Variant-level GxE analysis was implemented in GEM. LD clumping (Clump) and conditional and joint (COJO) analysis were performed using the combined GWAS of the main effect. Each dot indicates the performance metrics aggregating all three randomly selected loci, each having 50 simulation replicates under the corresponding simulation setting. Detailed statistics are provided in Supplementary Data 4. Calculations of power and empirical FDR for each method are described in Methods.With smaller effect size discrepancy (βe − βu), SharePro maintained a high power with a low FDR, whereas other methods had decreased power. For instance, in the simulation setting with three causal variants having a effect size discrepancy of 0.03 between the exposed and the unexposed category of equal sample size (βe  = 0.05; βu  = 0.0; Ne  = Nu = 25000 and KC = 3) (Fig. 3 and Supplementary Data4), SharePro achieved a power of 0.81 with an empirical FDR of 0.19 at the variant-level and a power of 0.98 with an empirical FDR of 0.02 at the effect group-level. After multiple testing correction, GEM had a power of 0.91 with an empirical FDR of 0.91 while both lead variants identified by COJO and LD-clumping achieved a power of 0.38 with an empirical FDR of 0.24 (Fig. 3 and Supplementary Data 4). This trend is similar with different number of causal variants (Supplementary Figs. 7–8 and Supplementary Data 4).SharePro differs from classical GxE methods in that it uses joint fine-mapping to reduce the multiple testing burden. We illustrate the importance of this difference with a simulation example (βe = 0.05; βu = 0.02; Ne = Nu = 25000, KC = 1) shown in Fig. 4. In this example, the association test in the exposed and the unexposed categories exhibited clear differences (Fig. 4A) and the GxE test in GEM accurately detected the simulated effect heterogeneity (Fig. 4B). However, the power for this test diminished after adjusting for multiple testing (Fig. 4B). In contrast, SharePro identified a candidate effect group through joint fine-mapping (Fig. 4C), thus avoiding testing on every variant, resulting in a successful GxE detection (Fig. 4D).Fig. 4: SharePro reduces multiple testing burden in GxE analysis.A Exposure stratified GWAS summary statistics in one simulation setting (βe  =  0.05; βu  = 0.02; Ne  = Nu  = 25, 000, KC  = 1). Each dot represents a variant and the color indicates its squared correlation with the causal variant. B GEM GxE p-values and Benjamini-Hochber (BH)-adjusted GEM GxE p-values. After multiple testing corrections, classical GxE analysis had insufficient power to detect effect heterogeneity. C Posterior inclusion probabilities in SharePro fine-mapping. Two variants in high LD were identified in an effect group. D The GxE p-value for the effect group identified in SharePro fine-mapping, represented by the top variant with the highest posterior inclusion probability in this effect group. All p-values are two-sided.SharePro identified genetic effects on lung function modulated by smoking statusBased on smoking status stratified GWAS for FFR, we investigated genetic effects on lung function modulated by smoking status. Across current smokers, past smokers and never smokers, we identified 375 effect groups associated with lung function (Supplementary Data 5). The genome-wide GxE analysis only detected genome-wide siginificant (p-value  < 5.0 × 10−8) effect heterogeneity mapped to the CHRNA3 gene comparing current smokers versus past smokers, current smokers versus never smokers, and past smokers versus never smokers (Fig. 5A). This GxSmoking interaction has been reported in literature15.Fig. 5: SharePro identifies genetic effects on lung function modulated by smoking status.A Manhattan plots demonstrating significance of genom-wide variant-level gene-by-smoking interactions (GxSmoking). The red dotted lines correspond to the genome-wide Bonferroni significance threshold. B Manhattan plots demonstrating significance of effect group-level GxSmoking. The red dotted lines correspond to the effect group-level Bonferroni significance threshold. Significant effect groups are labeled by their nearest genes. The unadjusted p-values are two-sided in (A) and (B). C Illustration of estimated effect sizes of effect groups that exhibited significant GxSmoking effects. The smoking status stratified GWAS included 27,409 current smokers, 100,285 past smokers and 155,951 never smokers. The effect groups are annotated using their nearest genes and top variants with the highest posterior inclusion probability in each effect group. The x-axis stands for estimated effect sizes and the colors represent different smoking statuses. The error bars represent 95% confidence intervals. FFR, the ratio of the forced expiratory volume in the first one second to the forced vital capacity of the lung measured by spirometry.With SharePro, we also detected an effect group mapped to the CHRNA3 gene (Fig. 5B), along with five additional effect groups with variants demonstrating different genetic effects on lung function in current smokers, past smokers and never smokers (Fig. 5C). In contrast, none of the lead variants detected by COJO and LD-clumping demonstrated significant effect heterogeneity after multiple testing correction (Supplementary Fig. 9).Interestingly, the estimated effect sizes of these effect groups across different smoking statuses suggest potential mechanisms of GxSmoking. The effect group mapped to the CHRNA3 gene demonstrated a large effect size in current smokers, a moderate effect size in past smokers and a minimal effect size in never smokers (Fig. 5C). The CHRNA3 gene is known to play a significant role in nicotine dependence and smoking behaviors16,17,18,19. Our results indicated that the genetic effect of CHRNA3 on lung function may be mediated by smoking, aligning with its biological function. Moreover, the effect groups mapped to the PRICKLE1 and the AMZ2 genes demonstrated association with FFR in current smokers but not in previous or never smokers, while the other effect groups exhibited more complex patterns of GxSmoking interaction (Fig. 5C). PRICKLE1 is an important regulator in the Wnt signaling pathway, which plays a crucial role in inflammation and immune response, as well as development and repair of lung tissues20,21,22. Although AMZ2 (archaelysin family metallopeptidase 2) has not been extensively studied, other human metalloproteases have been implicated in the regulation of airway inflammation23,24. Further investigations are necessary to confirm whether smoking-induced inflammation may disrupt such regulatory roles.We repeated the analyzes using inverse normal transformed FFR and obtained highly consistent results (Supplementary Fig. 10).SharePro improved characterization of sex differentiated genetic effects on fat distributionBased on sex stratified GWAS for WHRadjBMI, we investigated sex differentiated genetic effects on fat distribution and identified 334 effect groups associated with WHRadjBMI (Supplementary Data 6). Of these effect groups, 98 demonstrated effect group-level significance, including 47 demonstrating genome-wide significance (Supplementary Data 6).Effect groups demonstrated significant GxSex were mapped to genes that play important roles in adipose cell biology (Supplementary Data 6). For example, COBLL1 has been shown to play a crucial role in actin cytoskeleton remodelling during the differentiation of metabolically active and insulin-sensitive adipocytes25 and has been observed to exhibit sex dimorphic expression25. Recent research indicates RSPO3 can impact body fat distribution by regulating the expansion of gluteofemoral adipose tissue26 while estrogen is a key regulator of adipogenesis for gluteofemoral adipocytes27,28,29. Additionally, KLF14 is a transcription factor that regulates gene expression in adipose tissue30. Recent studies have demonstrated that KLF14 can increase adiposity and redistribute lipid storage in female mice, but not in male mice31. Moreover, SharePro identified an effect group with one variant rs1534696, mapped to the SNX10 gene, demonstrating significant GxSex in WHRadjBMI (p-value = 6.0 × 10−16; Supplementary Data 6). This variant has been shown to regulate diet-induced adipose expansion in female mice, but not in male mice32.We then assessed whether these sex differentiated genetic effects could be partially attributed to known factors of sex differences, including sex hormones and lipid metabolism. We observed that the top variants from effect groups demonstrating GxSex in WHRadjBMI tended to have sex differentiated effects on SHBG level (Spearman correlation = 0.42; p-value  = 1.5 × 10−15; Fig. 6A). In contrast, the correlation between sex differentiated genetic effects on WHRadjBMI and sex differentiated genetic effects on bioavailable testosterone was weak (Spearman correlation = 0.05; p-value = 0.35; Fig. 6B). Since SHBG regulates both testosterone and estrogen33, these results indicate that estrogen plays a more important role in sex differences of fat distribution. This aligns with the observation that most sex differentiated effect groups had a higher effect size in females (Supplementary Data 6) and supports the existing understanding that estrogen, rather than testosterone, promotes and maintains gluteofemoral fat distribution in females27,28,29. Furthermore, sex differentiated genetic effects on WHRadjBMI were more correlated with sex differentiated genetic effects on TG (Spearman correlation = 0.32; p-value = 3.5 × 10−9; Fig. 6C) than LDL (Spearman correlation = 0.24; p-value = 1.3 × 10−5Fig. 6D) and HDL (Spearman correlation = 0.16; p-value = 3.0 × 10−3Fig. 6E). Notably, effect groups mapped to hormone receptor coding genes with important roles in the control of feeding behaviors34,35,36, such as LHCGR, were effect group-level Bonferroni significant, which would not have been identified using the genome-wide significance threshold (Fig. 6C). After adjusting for levels of SHBG or lipids in the WHRadjBMI GWAS, the estimated GxSex interaction effects attenuated, though did not disappear (Supplementary Fig. 11), supporting that sex differentiated genetic effects on fat distribution may partially act through these molecules.Fig. 6: SharePro characterizes sex differentiated genetic effects on fat distribution.Effect size discrepancies in females and males for effect groups identified from sex stratified GWAS for waist-to-hip ratio adjusted for body mass index (WHRadjBMI) are plotted against (A) effect size discrepancies in females and males for sex hormone binding globulin (SHBG), (B) effect size discrepancies in females and males for bioavailable testosterone, (C) effect size discrepancies in females and males for triglycerides (TG), (D) effect size discrepancies in females and males for low-density lipoprotein (LDL) cholesterol, (E) effect size discrepancies in females and males for high-density lipoprotein (HDL) cholesterol, and (F) minor allele frequencies (MAF) discrepancies in females and males. The size of each dot corresponds to the unadjusted two-sided p-value of the effect group. Red dots represent genome-wide significant gene-by-sex interactions (GxSex), violet dots represent effect group-level Bonferroni significant GxSex, and gray dots represent other effect groups associated with WHRadjBMI. Effect groups reaching effect group-level significance are labeled by their nearest genes.Lastly, we did not observe sex differentiated minor allele frequencies in variants included in effect groups (Fig. 6F), indicating limited impact of sex differentiated participation bias37 in our study.
