Introduction

Polygenic scores (PGS) provide individualized genetic predictors of a phenotype by aggregating genetic effects across hundreds or thousands of loci, typically estimated from genome-wide association studies (GWAS). In recent years it has become increasingly apparent that the transferability of PGS performance across different cohorts is poor (1). Most analyses to date have focused on ancestry differences as the main driver of this lack of portability (2–4).

However, a growing body of evidence has demonstrated that PGS performance and effect estimates are influenced by differences in certain contexts i.e., environmental (classically termed “gene-environment” effects or interactions) or personal-level covariates – different phenotypes seem to be differently affected by these covariates, with adiposity traits such as body mass index (BMI) having substantial evidence for these effects (5–14). In one previous study, they showed that GWAS stratified by sample characteristics had better PGS performance in cohorts that matched the sample characteristics of the stratified GWAS, and that differences in heritability between the stratified cohorts partially explained this observation (13).

There are several gaps in current knowledge about these covariate-specific effects. Many analyses have assessed only a handful of these covariates due to the myriad of choices possible in typical large-scale biobanks. Little investigation has been done to systematically understand why certain covariates affect PGS performance, with such knowledge being useful to reduce the potential search for variables that impart context-specific effects. Furthermore, most studies investigating PGS-covariate interactions have been in European ancestry individuals; notably, comparing differences in PGS performance and prediction while controlling for differences in ancestry versus differences in context has not been assessed in previous studies. Moreover, covariate-specific effects are notorious for replicating poorly in human genetics studies, and previous studies of PGS-covariate interactions have been predominantly performed in the UK Biobank (UKBB) (15), where the majority of individuals are aged 40-69 (i.e., excluding young adults), are overall healthier than those from other e.g., hospital-based cohorts, and are predominantly European ancestry. Additionally, PGS performance is often assessed using linear models and in isolation of clinical covariates, which in practice would often be available.

Machine learning models can have increased performance over linear models and are capable of modeling complex relationships and interactions between variables, which may serve to increase predictive performance, especially given evidence for PGS-covariate specific effects. Finally, given evidence for context-specific effects, it should be possible to directly incorporate SNP-covariate interaction effects from a GWAS directly to improve prediction performance, instead of relying on post-hoc interactions from a typical PGS calculated from main GWAS effects.

Using genetic data with linked-phenotypic information from electronic health records, we estimated the effects of covariate stratification and interaction on performance and effect estimates of PGS for BMI (PGS_BMI) – a flowchart summarizing our analyses is presented in Figure 1. These analyses were done across four datasets (Supplemental Table 1): UK Biobank (UKBB), Penn Medicine BioBank (PMBB) (15), Electronic Medical Records and Genomics (eMERGE) network dataset (16), and Genetic Epidemiology Research on Adult Health and Aging (GERA). These datasets include participants from two ancestry groups (N=491,111 European ancestry (EUR), N=21,612 African ancestry (AFR)), and 62 covariates (25 present in multiple datasets) representing laboratory, survey, and biometric data types typically associated with cardiometabolic health and adiposity. After constructing PGS_BMI using out-of-sample multi-ancestry BMI GWAS, we assessed effects of covariate stratification on PGS_BMI R², the significance of PGS_BMI-covariate interaction terms and their increases to model R² over models only using main effects, as well as correlation of main effect, interaction effect, and R² differences. We then assessed ways to increase model performance through using machine learning models, and creating PGS_BMI using GxAge GWAS effects. This study addresses a plethora of open issues considering performance and effects of PGS on individuals from diverse backgrounds.

Figure 1.

Results

Effect of covariate stratification on PGS_BMI performance

We assessed 62 covariates for PGS_BMI R² differences (25 present, or suitable proxies, in multiple datasets (Supplemental Table 2) after stratifying on binary covariates and quintiles for continuous covariates. With UKBB EUR as discovery (N=376,729), 18 covariates had significant differences (Bonferroni p<.05/62) in R² among groups (Figure 2a), including age, sex, alcohol consumption, different physical activity measurements, Townsend deprivation index, different dietary measurements, lipids, blood pressure, and HbA1c, with 40 covariates having suggestive (p<.05) evidence of R² differences. From an original PGS_BMI R² of 0.076, R² increased to 0.094-0.088 for those in the bottom physical activity, alcohol intake, and high-density lipoprotein (HDL) cholesterol quintiles, and decreased to 0.067-0.049 for those in the top quintile, respectively, comparable to differences observed between ancestries (1). We note that the differences in R² due to alcohol intake and HDL were larger than those of any physical activity phenotype, despite physical activity having one of the oldest and most replicable evidence of interaction with genetic effects of BMI (17,18). Despite considerable published evidence suggesting covariate-specific genetic effects between BMI and smoking behaviors (6,8), we were only able to find suggestive evidence for R² differences when stratifying individuals across several smoking phenotypes (minimum p=0.016, for smoking pack years). R² differences due to educational attainment were also only suggestive (p=0.015), with published evidence on this association being conflicting (19–21).

Figure 2.

We replicated these analyses in three additional large-scale cohorts of European and African ancestry individuals (Figure 2b, Supplemental Table 3), as well as in African ancestry UKBB individuals. Among covariates with significant performance differences in the discovery analysis, we were able to replicate significant (p<.05) R² differences for age, HDL cholesterol, alcohol intake frequency, physical activity, and HbA1c, despite much smaller sample sizes. We again observed mostly insignificant differences across cohorts and ancestries when stratifying due to different smoking phenotypes and educational attainment. For each covariate and ancestry combination, we combined data across cohorts and conducted a linear regression weighted by sample size, regressing R² values on covariate values across groupings. Slopes of the regressions across cohorts had different signs between ancestries for the same covariate (triglyceride levels, HbA1c, diastolic blood pressure, and sex), although larger sample sizes may be needed to confirm these differences are statistically significant.

Several observations related to age-specific effects on PGS_BMI we considered noteworthy. First, in the weighted linear regression of all R² values across ancestries, expected R² for African ancestry individuals can become greater than that of European ancestry individuals among individuals within bottom and top age quintiles observed in these data. For instance, the predicted R² of 0.048 for 80 year-old European ancestry individuals would be lower than that of African ancestry individuals aged 24.7 and lower, indicating that differences in covariates can affect PGS_BMI performance more than differences due to ancestry. Second, we obtained these results despite the average age of GWAS individuals being 57.8, which should increase PGS_BMI R² for individuals closest to this age (13). This result suggests that PGS performance due to decreased heritability with age cannot be fully reconciled using GWAS from individuals of similar age being used to create PGS_BMI (as heritability is an upper bound on PGS performance). Finally, we observed that PGS_BMI R² increases as age decreases, consistent with published evidence suggesting that the heritability of BMI decreases with age (22,23).

PGS-covariate interaction effects

Next, we estimated difference in PGS effects due to interactions with covariates, by modeling interaction terms between PGS_BMI and the covariate for each covariate in our list (described in Methods). We implemented a correction for shared heritability between covariates of interest and outcome (which can inflate test statistics (24)) to better measure the environmental component of each covariate, and show that this correction successfully reduces significance of interaction estimates (Supplemental Figure 1). Again, using UKBB EUR as the discovery cohort, we observed 28 covariates with significant (Bonferroni p<.05/62) PGS-covariate interactions (Table 1), with 38 having suggestive (p<.05) evidence (Supplemental Table 4). We observed the largest effect of PGS-covariate interaction with alcohol drinking frequency (20.0% decrease in PGS effect per 1 standard deviation (SD) increase, p=2.62×10^-55), with large effects for different physical activity measures (9.4%-12.5% decrease/SD, minimum p=3.11×10^-66), HDL cholesterol (15.3% decrease/SD, p=1.71×10^-96) and total cholesterol (12.7% decrease/SD, p=1.64×10^-71). We observed significant interactions with diastolic blood pressure (10.8% increase/SD, p=6.06×10^{- 60}), but interactions with systolic blood pressure were much smaller (1.17% increase/SD, p=4.41×10^-3). Significant interactions with HbA1c (4.63% increase/SD, p=5.37×10^-14) and type 2 diabetes (27.2% PGS effect increase in cases, p=1.83×10^-7) were also observed. Other significant PGS-covariate interactions included lung function, age, sex, and LDL cholesterol – various dietary measurements also had significant interactions, albeit with smaller effects than other significant covariates. We were able to find significant interaction effects for smoking pack years (4.78% increase/SD, p=3.68×10^-7), but other smoking phenotypes had insignificant interaction effects after correcting for multiple tests (minimum p=2.7×10^-3); interactions with educational attainment were also insignificant (p= 4.54×10^-2).

Table 1.

We replicated these analyses across ancestries and the other non-UKBB EUR cohorts (Figure 3, Supplemental Table 4). For age and sex, which were available for all cohorts, interactions were significant (p<.05) and directionally consistent across cohorts and ancestries (except for GERA AFR which had small sample size (N=1,789)). We were able to test interactions with alcohol intake frequency and physical activity in GERA, and replicated significant and directionally consistent associations. We observed poor replication for LDL cholesterol, HbA1c, and smoking pack years, with insignificant and directionally inconsistent interaction effects across cohorts. Educational attainment was available in GERA, and interactions were once again insignificant. We observed significant and directionally consistent interaction effects for TG in eMERGE EUR and PMBB EUR, while the effect was inconsistent in UKBB EUR despite much larger sample size.

Figure 3.

However, despite significance of interaction terms, increases in model R² when including PGS-covariate interaction terms were small. For instance, the maximum increase among all covariates in UKBB EUR was only 0.0024 from a base R² of 0.1049 (2.1% relative increase), for alcoholic drinks per week. Across all cohorts and ancestries, the maximum increase in R² was only 0.0058 from a base R² of 0.09454 (6.1% relative increase), when adding a PGS-age interaction term for eMERGE EUR (p=5.40×10^-46) – this was also the largest relative increase among models with significant interaction terms. This result suggests that, while interaction effects can significantly modify PGS_BMI effect, their overall impact on model performance is relatively small, despite large differences in R² when stratifying by covariates.

Correlations between R² differences, interaction effects, and main effects

We next investigated the relationship between interaction effects, maximum R² differences across quintiles, and main effects of covariates on BMI. We first estimated main effects of each covariate on BMI (Methods, Supplemental Table 5), and then calculated the correlation weighted by sample size between main effects, maximum PGS_BMI R² across quintiles, and PGS-covariate interaction effects (Figure 4) across all cohorts and ancestries – GERA data were excluded from these analyses due to slightly different phenotype definition (Supplemental Table 6), as were binary variables. Interaction effects and maximum R² differences had a 0.80 correlation (p=2.1×10^-27), indicating that variables with larger interaction effects also had larger effects on PGS_BMI performance across quintiles, and that covariates that increase PGS_BMI effect also have the largest effect on PGS_BMI performance i.e., individuals most at risk for obesity will have both disproportionately larger PGS_BMI effect and R². Main effects and maximum R² differences had a 0.56 correlation (p=1.3×10^-11), while main effects and interaction effects had a 0.58 correlation (p=7.6×10^-12) again suggesting that PGS_BMI are more predictive in individuals with higher values of BMI-associated covariates, although less predictive than estimating the interaction effects themselves directly. However, this result demonstrates that covariates that influence both PGS_BMI effect and performance can be predicted just using main effects of covariates, which are often known for certain phenotypes and easier to calculate, as genetic data and PGS construction would not be required.

Figure 4.

Increase in PGS effect for increasing percentiles of BMI itself, and its relation to R² differences when stratifying by covariates Given large and replicable correlations between main effects, interaction effects, and maximum R² differences for individual covariates, it seemed these differences may be due to the differences in BMI itself, rather than any individual or combination of covariates. To assess this, we used quantile regression to evaluate the effect of PGS_BMI on BMI at different deciles of BMI itself. We observed that the effect of PGS_BMI consistently increases from lower deciles to higher deciles across all cohorts and ancestries (Figure 5) – for instance, in European ancestry UKBB individuals, the effect of PGS_BMI (in units of log(BMI)) when predicting the bottom decile of log(BMI) was 0.716 (95% CI: 0.701-.732), and increased to 1.31 (95% CI: 1.29-1.33) in the top decile. Across all cohorts and ancestries, the effect of PGS_BMI between lowest and highest effect decile ranged from 1.43-2.06 times larger, with all cohorts and ancestries having non-overlapping 95% confidence intervals between their effects (except for African ancestry eMERGE individuals, which had much smaller sample size).

Figure 5.

While this analysis showed that the effect of PGS_BMI increases as BMI itself increases, which may help explain significant interaction effects between PGS_BMI and different covariates, it does not directly explain differences in R² when stratifying by different covariates – we describe several points that help explain this result and suggest they may actually be closely related. Essentially, as the magnitude of the slope of a regression line increases while the mean squared residual does not increase, model R² will increase – we demonstrate this using simulated data (Supplemental Figure 2). As the magnitude of the regression line’s slope decreases, the regression line becomes a comparatively worse predictor compared to just using the mean, which decreases R² despite the mean error being the same across models. To demonstrate this in real data, we compared simple univariable models of log(BMI) ∼ PGS_BMI (in units of log(BMI)) between the bottom and top age quintiles in the European ancestry UKBB (Supplemental Figure 3). As shown in previous sections, R² and PGS_BMI beta are higher in younger individuals (R² = 0.088 versus R² = 0.066, and beta=1.12 and 0.87, respectively), which seem to be a direct consequence of one another, as the mean squared error in younger individuals is actually higher (0.027 versus 0.022, respectively). This description suggests that the use of R² as the sole performance metric for evaluation of PGS may not always be appropriate, despite its overwhelming usage. Furthermore, this explanation helps explain the seemingly paradoxical results of significant interaction terms yet small increases in overall model R², and comparably much larger differences in R² in the stratified analyses.

S Figure 2.

S Figure 3.

Effects of machine learning approaches on predictive performance

Given evidence of PGS-covariate dependence, we aimed to assess increases in R² when using machine learning models (neural networks), which can inherently model interactions and other nonlinearities, over linear models even with interaction terms. We first included age and sex as the only covariates (along with genotype PCs and PGS_BMI), as age and sex were present in all datasets and had significant and replicable evidence for PGS-dependence across our analyses. Three models were assessed – L1-regularized (i.e., LASSO) linear regression without any interaction terms, LASSO including a PGS-age and PGS-sex interaction term, and neural networks (without interaction terms). When comparing neural networks to LASSO with interaction terms, the relative 10-fold cross-validated R² increased up to 67% (mean 23%) across cohorts and ancestries (Figure 6, Supplemental Table 7). The inclusion of interaction terms increased cross-validated R² up to 12% (mean 3.9%) when comparing LASSO including interaction terms to LASSO with main effects only.

Figure 6.

We then modeled all available covariates and their interactions with PGS for each cohort and did similar comparisons. Cross-validated R² increased by up to 17.6% (mean 9.5%) when using neural networks versus LASSO with interaction terms, and up to 7.0% (mean 2.0%) when comparing LASSO with interaction terms to LASSO with main effects only. Increases in model performance using neural networks were smaller in UKBB, perhaps due to the age range being smaller than other cohorts (all participants aged 39-73). This result suggests that additional variance explained through non-linear effects with age and sex are explained by other variables present in the remainder of the datasets. Our findings show machine learning methods can improve model R² that include PGS_BMI as variables beyond including interaction terms in linear models, even when variable selection is performed using LASSO, demonstrating that model performance can be increased beyond modeling nonlinearities through linear interaction terms and a feature selection procedure.

Calculating PGS directly from GxAge GWAS effects

Previous studies (13) have created PGS using GWAS stratified by different personal-level covariates, but for practical purposes this leads to a large loss of power as the full size of the GWAS is not utilized for each strata and continuous variables are forced into bins. We developed a novel strategy where PGS are instead created from a full-size GWAS that includes SNP-covariate interaction terms (Methods). We focused on age interactions, given their large and replicable effects based on our results – similar to a previous study (13), we conducted these analyses in the European UKBB. We used a 60% random split of study individuals to conduct three sets of PGS using GWAS of the following designs: main effects only, main effects also with a SNP-age interaction term, and main effects but stratified into quartiles by age. 20% of the remaining data were used for training and the final 20% were held-out as a test set. The best performing PGS created from SNP-age interaction terms (PGS_GxAge) increased test R² to 0.0771 (95% bootstrap CI: 0.0770-0.0772) from 0.0715 (95% bootstrap CI: 0.0714-0.0716), a 7.8% relative increase compared to the best performing main effect PGS (Figure 7, Supplemental Table 8 – age-stratified PGS had much lower performance than both other strategies (unsurprising given reasons previously mentioned). Including a PGS_GxAge-age interaction term only marginally increased R² (0.0001 increase), with similar increases for the other two sets of PGS, further demonstrating that post-hoc modeling of interactions cannot reconcile performance gained through directly incorporating interaction effects from the original GWAS. The strategy of creating PGS directly from full-sized SNP-covariate interactions is potentially quite useful as it increases PGS performance without the need for additional data – there are almost certainly a variety of points of improvement (described more in Discussion), but we consider their investigation outside the scope of this study.

Figure 7.

Discussion

We uncovered replicable effects of covariates across four large-scale cohorts of diverse ancestry, on both performance and effects of PGS_BMI. When stratifying by quintiles of different covariates, certain covariates had significant and replicable evidence for differences in PGS_BMI R², with R² being nearly double between top and bottom performing quintiles for covariates with the largest differences. When testing PGS-covariate interaction effects, we also found covariates with significant interaction effects, where, for the largest effect covariates, each standard deviation change affected PGS_BMI effect by nearly 20%. Across analyses, we found age and sex had the most replicable interaction effects, with levels of serum cholesterol, physical activity, and alcohol consumption having the largest effects across cohorts. Interaction effects and R² differences were strongly correlated, with main effects also being correlated with interaction effects and R² differences, suggesting that covariates with the largest interaction effects also contribute to the largest R² differences, with simple main effects also being predictive of expected differences in R² and interaction effects. Relatedly, we observed the effect of PGS_BMI increases as BMI itself increases, and reason that differences in R² when stratifying by covariates are largely a consequence of difference in PGS_BMI effects. Next, we employed machine learning methods for prediction of BMI with models that include PGS_BMI and demonstrate that these methods outperform regularized linear regression models that include interaction effects. Finally, we employed a novel strategy that directly incorporates SNP interaction effects into PGS construction, and demonstrate that this strategy improves PGS performance when modeling SNP-age interactions compared to PGS created only from main effects.

These observations are relevant to current research and clinical use of PGS, as individuals above a percentile cutoff are designated high-risk (40), implying that individuals most at-risk for obesity have both disproportionately higher predicted BMI and increased BMI prediction performance compared to the general population. More broadly, these results may likely extend to single variant effects instead of those aggregated into a PGS, which may inform the cause of previous GxE discoveries – for instance, variants near FTO that interact with physical activity discovered through GWAS of BMI are robust and well-documented. However, individuals engaging in physical activity will generally have lower BMI than those that are sedentary, and these results suggest it may not be the difference in physical activity that’s driving the interaction, rather the difference in BMI itself. This concept may also apply to other traits – for instance, sex-specific analyses are commonly performed, and variants with differing effects between males and females GWAS may largely be explained by phenotypic differences, rather than any combination of biological or lifestyle differences.

Future work may include replicating these analyses across additional traits, and trying to understand why these differences occur, as well as further exploring machine learning and deep learning methods on other phenotypes to determine if this trend of inclusion of PGS along with covariate interaction effects outperforms linear models for risk prediction. Additionally, inclusion of a PGS for the covariate to better measure its environmental effect is potentially worth exploring further, and should improve in the future as PGS performance continues to increase. A slight limitation of this method in our study is that for the UKBB analyses the GWAS used for PGS construction were also from UKBB thus not out-of-sample, although many of the covariates only have GWAS available through UKBB individuals. Furthermore, a variety of improvements are likely possible when creating PGS directly from SNP-covariate interaction terms. First, we used the same SNPs that were selected by pruning and thresholding based on their main effect p-values, but selection of SNPs based on their interaction p-values should also be possible and would likely improve performance. Additionally, performance of pruning and thresholding-based strategies have largely been overtaken by methods that first adjust all SNP effects for LD and don’t require exclusion of SNPs, and a method that could do a similar adjustment for interaction effects would likely outperform most current methods for traits with significant context-specific effects. Next, incorporating additional SNP-covariate interactions (e.g., SNP-sex) would also likely further improve prediction performance, although any SNP selection/adjustment procedures may be further complicated by additional interaction terms.

Finally, if SNP effects do truly differ according to differences in phenotype, then SNP effects would differ depending on the alleles one has, implying epistatic interactions are occurring at these SNPs.

While difference in phenotype itself may be able to explain difference in genetic effects, it still may be that specific environmental or lifestyle characteristics are driving the differences. We propose several ideas about why BMI-associated covariates have larger interaction effects and impact on R² for PGS_BMI. First, age may be a proxy for accumulated gene-environment interactions as younger individuals have less exposure to environmental influences on weight compared to older individuals; therefore, one may expect that in younger individuals their phenotype could be better explained by genetics compared to environment. Second, PGS may more readily explain high phenotype values especially for positively skewed phenotypes, as large effect variants (e.g., associated with very high weight or height (25)) may be more responsible for extremely high phenotypic values. For example, the distribution of BMI is often positively skewed, and effects in trait-increasing alleles may have a larger potential to explain trait variation compared to trait-reducing variants. This explanation would likely be better suited to positively skewed traits and is not fully satisfactory as first log-transforming or rank-normal transforming the phenotype, as was done in this study, may invalidate this explanation.

PGS is a promising technique to stratify individuals for their risk of common, complex disease. To achieve more accurate predictions as well as promote equity, further research is required regarding PGS methods and assessments. This research provides firm evidence supporting the context-specific nature of PGS and the impact of nonlinear covariate effects for improving polygenic prediction of BMI, promoting equitable use of PGS across ancestries and cohorts.

Methods

Study datasets

Individual inclusion criteria and sample sizes per cohort are described below – additional information is available in Supplemental Table 1.

UKBB

Individual-level quality-control and filtering have been described elsewhere (26) for European ancestry individuals. Briefly, individuals were split by ancestry according to both genetically inferred ancestry and self-reported ethnicity (15). Individuals with low genotyping quality and sex mismatch were removed, only unrelated individuals (pi-hat < 0.250) were retained, and variants were filtered at INFO > 0.30 and minor allele frequency (MAF) > 0.01. For African ancestry, individuals were first selected based on self-reported ethnicity “Black or Black British”, “Caribbean”, “African”, or “Any other Black background”. Individuals that were low quality i.e., “Outliers for heterozygosity or missing rate”, and that were Caucasian from “Genetic ethnic grouping” were removed. Of these individuals, those that were +/- 6 standard deviations from the mean of the first 5 genetic principal components provided by UKBB were excluded.

Finally, only unrelated individuals were retained up to the second degree using plink2 (27) “-king-cutoff 0.125”. After QC and consideration of phenotype, a total of 7,046 individuals in the UKBB AFR data who also had BMI available were used for downstream analyses. In total, 383,775 individuals were used for analysis (N_EUR=376,729, N_AFR=7,046).

eMERGE

Ancestry and relatedness inference have been described elsewhere (15). Individuals were split into European and African ancestry cohorts, and related individuals were removed (pi-hat > 0.250) such that all were unrelated. 35,064 individuals (N_EUR=31,961, N_AFR=3,103) were used for analysis.

GERA

Ancestry inference has been described elsewhere (28), and study individuals were divided into European and African ancestry cohorts. Related individuals were removed using plink2 “-king-cutoff 0.125”. 57,838 individuals (N_EUR=56,049, N_AFR=1,789) were used for analysis.

PMBB

Ancestry inference and relatedness inference have been described elsewhere (29). Individuals were split into European and African ancestry cohorts, and related individuals were removed at pi-hat > 0.250. 36,046 individuals (N_EUR=26,372, N_AFR=9,674) were used for analysis.

Choice of covariates

A total of 62 covariates were included in the analyses, 25 of which were present (or similar proxies) in multiple datasets. These covariates were selected based on relevance to cardiometabolic health and obesity, and previous evidence of context-specific effects with BMI (5,6,8,30–32). For UKBB, phenotype values were used from the collection that was closest to recruitment, and for PMBB the median values were used – for GERA and eMERGE only one value was available. Additional details on covariate constructions, transformations, filtering, and cohort availability are in Supplemental Table 2.

PGS construction

PGS for BMI (PGS_BMI) were constructed using PRS-CSx (33), using GWAS summary statistics for individuals of European (34), African (35), and East Asian (36) ancestry that were out-of-sample of study participants. A set of 1.29 million HapMap3 SNPs provided by PRS-CSx was used for PGS calculation, which are generally well-imputed and variable frequency across global populations. Default settings for PRS-CSx (downloaded November 22, 2021) were used, which have been demonstrated to perform well for highly polygenic traits such as BMI (list of parameters in Supplemental Table 9). The final PGS_BMI per ancestry and cohort was calculated by regressing log(BMI) on the PGS_BMI per ancestry without covariates – the combined, predicted value was used as a single PGS_BMI in downstream analyses.

For GERA, BMI was not transformed, as it was already binned into a categorical variable with 5 levels (18≤, 19-25, 26-29, 30-39, >40). Additionally, for GERA the uncombined ancestry-specific PGS_BMI were used in the final models, as it had higher R² than using the combined PGS_BMI (data not shown).

PGS_BMI performance after covariate stratification

Analyses were performed separately for each cohort and ancestry. For each covariate, individuals were binned by binary covariates or quintiles for continuous covariates. Incremental PGS_BMI R² was calculated by taking the difference in R² between:

We performed 5,000 bootstrap replications to obtain a bootstrapped distribution of R². P-values for differences in R² were calculated between groups by calculating the proportion of overlap between two normal distributions of the R² value using the standard deviations of the bootstrap distributions. Again for GERA, BMI was not transformed.

PGS_BMI interaction modelling

Evidence for interaction with each covariate with the PGS_BMI was evaluated using linear regression. It has been reported that the inclusion of covariates that are genetically correlated with the outcome can inflate test statistic estimates (24,37,38). To assuage these concerns, we introduced a novel correction, where we first calculated a PGS for the covariate (PGS_Covariate) and included it in the model, as well as an interaction term between PGS_BMI and PGS_Covariate. The PGS_Covariate terms were calculated using the European ancestry Neale Lab summary statistics (URLs) and PRS-CS (39). To standardize effect sizes across analyses, PGS_BMI and Covariate were first converted to mean zero and standard deviation of 1 (binary covariates were not standardized). We demonstrate inclusion of PGS_Covariate terms successfully reduced significance of the PGS_BMI*Covariate term (Supplemental Figure 1). The final model used to evaluate PGS_BMI and Covariate interactions was:

S Figure 1.

We report the effect estimates of the PGS_BMI*Covariate term, and differences in model R² with and without the PGS_BMI*Covariate term. Again for GERA, BMI was not transformed.

Correlation between R² differences, interaction effects, and main effects

We estimated main effects of each covariate on BMI with the following model:

Note that we ran new models with main effects only, instead of using the main effect from the interaction models (as the main effects in the interaction models depend on the interaction terms, and main effects used to create interaction terms are sensitive to centering of variables despite the scale invariance of linear regression itself (40)). We then estimated the correlation between main effects, interaction effects, and maximum R² differences across all cohorts and ancestries weighting by sample size, analyzing quantitative and binary variables separately.

Quantile regression to measure PGS effect across percentiles of BMI

The effect of PGS_BMI on BMI at different deciles of BMI was assessed using quantile regression. Tau – the parameter that sets which percentile to be predicted – was set to .10, .20, …, .90.

Models included age, sex, and the top 5 genetic PCs as covariates. Analyses were stratified by ancestry and cohort, and BMI was first log transformed. GERA was excluded from these analyses, as a portion of the models failed to run (as BMI values from GERA were already binned, some deciles all had the same BMI value – additionally, difference in effects between bins would be harder to evaluate as BMI within each decile would be more homogeneous).

Machine learning models

UKBB EUR and GERA EUR models were restricted to 30,000 random individuals, for computational reasons – BMI distributions did not differ from the full-sized datasets (Kolmogorov-Smirnov p-value of 0.29 and 0.57, respectively). PGS_BMI and top five genetic principal components were included as features in all models. Two sets of models were evaluated for each cohort and ancestry: including age and sex as features, and including all available covariates in each cohort as features. Interactions terms between PGS and each covariate were included for models using interaction terms. ‘Ever Smoker’ status was used in favor of ‘Never’ vs ‘Current smoking’ status (if present), as individuals with ‘Never’ vs ‘Current’ status are a subset of those with ‘Ever Smoker’ status. UKBB AFR with all covariates was excluded due to small sample size (N=53).

Neural networks were used as the model of choice, given their inherent ability model interactions and other nonlinear dependencies. Prior to modeling, all features were scaled to be between 0 and 1. We used average 10-fold cross validation R² to evaluate model performance. Separate models were trained using untransformed and log(BMI). L1-regularized linear regression was used with 18 values of lambda (1.0, 5.0×10^-1, 2.0×10^-1, 1.0×10^-1, 5.0×10^-2, 2.0×10^{- 2}, …, 1.0×10^-5, 5.0×10^-6, 2.0×10^-6). Models were trained without inclusion of interaction terms (which neural networks can implicitly model), using 1,000 iterations of random search with the following hyperparameter ranges: size of hidden layers [10, 200], learning rate [0.01, 0.0001], type of learning rate [constant, inverse scaling], power t [0.4, 0.6], momentum [0.80, 1.0], batch size [32, 256], and number of hidden layers [1, 2].

GxAge PGS_BMI creation and assessment

Analyses were conducted in the European UKBB (N=376,629), as was done in a study on a similar topic (13). Three sets of analyses were performed, using GWAS conducted in a 60% random split of individuals using the following models (BMI was rank-normal transformed before GWAS):

1. BMI ∼ SNP + Age + Sex + PCs_1-5

2. BMI ∼ SNP + Age*SNP + Age + Sex + PCs_1-5

3. Using the model in 1) but stratified into quartiles by age. BMI was rank-normal transformed within each quartile.

Using each set of GWAS, PGS was first calculated in a 20% randomly selected training set of the dataset using pruning and thresholding using 10 p-value thresholds (0.50, 0.10, …, 5.0×10^-5, 1.0×10^-5) and remaining settings as default in Plink 1.9. For 2), GxAge PGS_BMI were calculated using SNPs clumped by their main effect p-values from 1), and additionally incorporating the GxAge interaction effects per SNP. In other words, instead of typical PGS construction as:

PGS_i = β₁k₁ + β₂k₂ + … + β_nk_n

For an individual i’s PGS calculation, with main SNP effect β, and n SNPs, PGS incorporating GxAge effects (PGS_GxAge) was calculated as:

Where β_GxAge is the GxAge effect for each SNP n and Age_i is the age for individual i.

For each of the three analyses, the parameters and models resulting in the best performing PGS (highest incremental R², using same main effect covariates as in the three GWAS) from the training set were evaluated in the remaining 20% of the study individuals. For 3), models were first trained within each quartile separately. To maintain sense of scale across quartiles (after rank-normal transformation), R² between all predicted values and true values was calculated together. For R² confidence intervals, the training set was bootstrapped and evaluated on the test set 5,000 times.

URLs

Neale Lab UKBB summary statistics: http://www.nealelab.is/uk-biobank

Regeneron Genetics Center Banner Author List and Contribution Statements

RGC Management and Leadership Team

Goncalo Abecasis, PhD, Aris Baras, M.D., Michael Cantor, M.D., Giovanni Coppola, M.D., Andrew Deubler, Aris Economides, Ph.D., Luca A. Lotta, M.D., Ph.D., John D. Overton, Ph.D., Jeffrey G. Reid, Ph.D., Katherine Siminovitch, M.D., Alan Shuldiner, M.D.

Sequencing and Lab Operations

Christina Beechert, Caitlin Forsythe, M.S., Erin D. Fuller, Zhenhua Gu, M.S., Michael Lattari, Alexander Lopez, M.S., John D. Overton, Ph.D., Maria Sotiropoulos Padilla, M.S., Manasi Pradhan, M.S., Kia Manoochehri, B.S., Thomas D. Schleicher, M.S., Louis Widom, Sarah E. Wolf, M.S., Ricardo H. Ulloa, B.S.

Clinical Informatics

Amelia Averitt, Ph.D., Nilanjana Banerjee, Ph.D., Michael Cantor, M.D., Dadong Li, Ph.D., Sameer Malhotra, M.D., Deepika Sharma, MHI, Jeffrey Staples, Ph.D.

Genome Informatics

Xiaodong Bai, Ph.D., Suganthi Balasubramanian, Ph.D., Suying Bao, Ph.D., Boris Boutkov, Ph.D., Siying Chen, Ph.D., Gisu Eom, B.S., Lukas Habegger, Ph.D., Alicia Hawes, B.S., Shareef Khalid, Olga Krasheninina, M.S., Rouel Lanche, B.S., Adam J. Mansfield, B.A., Evan K. Maxwell, Ph.D., George Mitra, B.A., Mona Nafde, M.S., Sean O’Keeffe, Ph.D., Max Orelus, B.B.A., Razvan Panea, Ph.D., Tommy Polanco, B.A., Ayesha Rasool, M.S., Jeffrey G. Reid, Ph.D., William Salerno, Ph.D., Jeffrey C. Staples, Ph.D., Kathie Sun, Ph.D.

Analytical Genomics and Data Science

Goncalo Abecasis, D.Phil., Joshua Backman, Ph.D., Amy Damask, Ph.D., Lee Dobbyn, Ph.D., Manuel Allen Revez Ferreira, Ph.D., Arkopravo Ghosh, M.S., Christopher Gillies, Ph.D., Lauren Gurski, B.S., Eric Jorgenson, Ph.D., Hyun Min Kang, Ph.D., Michael Kessler, Ph.D., Jack Kosmicki, Ph.D., Alexander Li, Ph.D., Nan Lin, Ph.D., Daren Liu, M.S., Adam Locke, Ph.D., Jonathan Marchini, Ph.D., Anthony Marcketta, M.S., Joelle Mbatchou, Ph.D., Arden Moscati, Ph.D., Charles Paulding, Ph.D., Carlo Sidore, Ph.D., Eli Stahl, Ph.D., Kyoko Watanabe, Ph.D., Bin Ye, Ph.D., Blair Zhang, Ph.D., Andrey Ziyatdinov, Ph.D.

Therapeutic Area Genetics

Ariane Ayer, B.S., Aysegul Guvenek, Ph.D., George Hindy, Ph.D., Giovanni Coppola, M.D., Jan Freudenberg, M.D., Jonas Bovijn M.D., Katherine Siminovitch, M.D., Kavita Praveen, Ph.D., Luca A. Lotta, M.D., Manav Kapoor, Ph.D., Mary Haas, Ph.D., Moeen Riaz, Ph.D., Niek Verweij, Ph.D., Olukayode Sosina, Ph.D., Parsa Akbari, Ph.D., Priyanka Nakka, Ph.D., Sahar Gelfman, Ph.D., Sujit Gokhale, B.E., Tanima De, Ph.D., Veera Rajagopal, Ph.D., Alan Shuldiner, M.D., Bin Ye, Ph.D., Gannie Tzoneva, Ph.D., Juan Rodriguez-Flores, Ph.D.

Research Program Management & Strategic Initiatives

Esteban Chen, M.S., Marcus B. Jones, Ph.D., Michelle G. LeBlanc, Ph.D., Jason Mighty, Ph.D., Lyndon J. Mitnaul, Ph.D., Nirupama Nishtala, Ph.D., Nadia Rana, Ph.D., Jaimee Hernandez

Data Availability

UK Biobank data was accessed under project #32133. eMERGE data is available at dbGaP in phs001584.v2.p2. GERA data is available at dbGaP in phs000674.v3.p3. Summary statistics for PMBB are available from authors upon request.

Acknowledgements

We acknowledge David Crosslin for helping clean the eMERGE data.

Funding

M.D.R., S.D., and D.H. are funded by AI077505 and HL169458. W.K.C is funded by grant NIDDK DK52431. J.F.P. is funded by grant U01 HG011166. C.W. is funded by U01 HG008680. This phase of the eMERGE Network was initiated and funded by the NHGRI through the following grants: U01HG008657 (Group Health Cooperative/University of Washington); U01HG008685 (Brigham and Women’s Hospital); U01HG008672 (Vanderbilt University Medical Center); U01HG008666 (Cincinnati Children’s Hospital Medical Center); U01HG006379 (Mayo Clinic); U01HG008679 (Geisinger Clinic); U01HG008680 (Columbia University Health Sciences); U01HG008684 (Children’s Hospital of Philadelphia); U01HG008673 (Northwestern University); U01HG008701 (Vanderbilt University Medical Center serving as the Coordinating Center); U01HG008676 (Partners Healthcare/Broad Institute); U01HG008664 (Baylor College of Medicine); and U54MD007593 (Meharry Medical College).