Introduction

Polygenic scores (PGS) provide individualized genetic predictors of a phenotype by aggregating genetic effects across hundreds or thousands of loci, typically 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 transferability (24). However, a growing body of evidence has demonstrated that PGS performance and effect estimates are influenced by differences in certain environmental (classically termed “gene-environment” effects or interactions) or personal-level covariates – different phenotypes seem to be differently affected by these covariates, with most evidence being for adiposity traits such as body mass index (BMI) (514).

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; unfortunately, 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 having poor replicability in human genetics studies, and previous PGS-covariate interaction studies 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.

In this study we aimed to address the previously mentioned concerns. 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 (PGSBMI) – Figure 1 is a flowchart of the project. These analyses were done across four datasets: 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 (including blood lipids, socioeconomic status, blood pressure, diet, physical activity, alcohol intake, smoking, and lung function, among others). We calculated the difference in PGSBMI R2 among individuals stratified by different quintiles of these variables (different groups for binary variables) and assessed whether a variable’s association with BMI could explain its effect on PGSBMI performance across groups. We also assessed the significance of PGSBMI-covariate interaction terms, and their increases to model R2 over models only using main effects, as well as correlation of main effect, interaction effect, and R2 differences. Replication and comparisons across all datasets and ancestries were conducted. We then used machine learning models and evaluated their increase in performance over linear models, even those that included regularization and interaction terms. Finally, we conducted BMI GWAS using GxAge interaction terms, created PGSBMI directly using these effects, and assessed performance gains over using typical main effect GWAS. This study addresses a plethora of open issues considering performance and effects of PGS on individuals from diverse backgrounds.

A flowchart of the project.

Results

Effect of covariate stratification on PGSBMI performance

We assessed 62 covariates for R2 differences (25 present, or suitable proxies, in multiple datasets) 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 R2 among groups (Figure 2a), including age, sex, alcohol consumption, different physical activity measurements, Townsend deprivation index, different dietary measurements, lipids, blood pressure, and A1c, with 40 covariates having suggestive (p<.05) evidence of R2 differences. From an original PGSBMI R2 of .076, R2 increased to .094, .093, and .088 for those in the bottom physical activity, alcohol intake, and high-density lipoprotein (HDL) cholesterol quintiles, and decreased to .067, .57, and .049 for those in the top quintile (71%, 61%, and 56% relative R2), respectively, comparable to differences due to ancestry (1). We note that the differences in R2 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 (33,34). 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 R2 differences when stratifying individuals across several smoking phenotypes (minimum p=0.016, for smoking pack years). R2 differences due to educational attainment were also only suggestive (p=0.015), of which published evidence is conflicting (35– 37).

PGS R2 stratified by quintiles for quantitative variables and by binary variables. a) Continuous covariates with significant (p < 7.7×10−4) R2 differences across quintiles in UKBB EUR. Pork and processed meat consumption per week were excluded from this plot in favor of pork and processed meat intake. b) Covariates with significant differences that were available in multiple cohorts, note actual values instead of percentiles plotted on x-axis (except percentiles were used for physical activity, alcohol intake frequency, and socioeconomic status, which had slightly differing phenotype definitions across cohorts). Townsend index and income were used as variables for socioeconomic status UKBB and GERA, respectively. Note that the sign for Townsend index was reversed, since increasing Townsend index is lower socioeconomic status, while increasing income is higher socioeconomic status. Abbreviations: physical activity (PA), International Physical Activity Questionnaire IPAQ).

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. For covariates with significant performance differences in the discovery analysis, we were able to replicate significant (p<.05) R2 differences for age, HDL cholesterol, alcohol intake frequency, physical activity, and HbA1c, despite much smaller sample sizes. We again observed largely insignificant differences across cohorts and ancestries when stratifying due to different smoking phenotypes and educational attainment. Combining across cohorts, for each covariate and ancestry we regressed R2 values across groupings on covariate values in a weighted linear regression, weighting by sample size. Slope of a weighted linear regression across cohorts was directionally inconsistent 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.

We noted several observations related to age-specific effects on PGSBMI. First, in the weighted linear regression of all R2 across ancestries, expected R2 for African ancestry individuals overcomes that of European ancestry individuals among individuals within bottom and top age quintiles observed in these data. For instance, the predicted R2 of .048 for 80 year-old European ancestry individuals would be lower than that of African ancestry individuals aged 24.7 and lower, indicating that certain differences in covariates can overcome differences in PGSBMI performance due to ancestry (with age being one of the most common variables used for assessing and determining health). Second, we obtained these results despite the average age of GWAS individuals being 57.8, which should increase PGSBMI R2 of individuals closest to this age (13). This result suggests that in the case of these analyses, PGS performance due to decreased heritability with age cannot be reconciled by having GWAS individuals of similar age being used to create the PGSBMI. Finally, we observed that PGSBMI R2 increases as age decreases, consistent with published evidence suggesting that the heritability of BMI decreases with age (38,39).

PGS-covariate interaction effects

Next, we estimated difference in PGS effects due to interactions with covariates, by modeling interaction terms between PGSBMI for each covariate (Methods). We implemented a correction for shared heritability between covariates of interest and outcome (which can inflate test statistics (28)) 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 discovery, 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 interactions after multiple testing correction (minimum p=2.7×10−3) – interactions with educational attainment were also insignificant (p= 4.54×10−2).

Model descriptive statistics on 28 of 62 covariates, which have significant (p<.05/62) PGS-covariate interaction terms, in UKBB EUR. The third column is the percentage change in PGS effect per unit change (standard deviations for continuous variables, binary variables encoded as 0 or 1) in covariate. The fifth column is the increase in model R2 with a PGS-covariate interaction term versus a main effects only model. Abbreviations: blood pressure (BP), physical activity (PA), forced vital capacity (FVC), forced expiratory volume in 1-second (FEV1), International Physical Activity Questionnaire (IPAQ).

We replicated these analyses across ancestries and the other cohorts in this study (Figure 3, Supplemental Table 4). For age and sex, which were available for all cohorts, interactions were significant (p<.05) and directionally consistent across all cohorts and ancestries (except for GERA AFR both effects were reversed, albeit with small sample size (N=1,789)). We were able to test interactions with alcohol intake frequency and physical activity in GERA, and we were able to replicate significant and directionally consistent association. 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.

Relative percentage changes in PGS effect per unit change in covariate, for covariates that significantly changed PGS effect (i.e., significant interaction beta at Bonferroni p < 7.7×10−4– denoted by asterisks) and were present in multiple cohorts and ancestries. Same covariate groupings and transformations were performed as with Figure 1.

However, despite significance of interaction terms, increases in model R2 when including PGS-covariate interaction terms were small. For instance, the maximum increase among all covariates in UKBB EUR was only .0024 from a base R2 of .1049 (2.1% relative increase), for alcohol drink frequency per week. Across all cohorts and ancestries, the maximum increase in R2 was only .0058 from a base R2 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 the effect of PGSBMI, their overall impact on model performance is relatively small, despite the observation that stratifying by covariates produces large differences in R2.

Correlations between R2 differences, interaction effects, and main effects

We next investigated the relationship between interaction effects, maximum R2 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 PGSBMI R2 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 R2 differences had a .80 correlation (p=2.1×10−27), indicating that variables with the largest interaction effects also had the largest effect on PGSBMI performance across quintiles, and that covariates that increase PGSBMI effect also have the largest effect on PGSBMI performance i.e., individuals most at risk for obesity will have both disproportionately larger PGSBMI effect and R2. Main effects and maximum R2 differences had a .56 correlation (p=1.3×10−11), while main effects and interaction effects had a .58 correlation (7.6×10−12) again suggesting that PGSBMI are more predictive in individuals with higher values of BMI-associated covariates, although predictive to a lower degree than estimating the interaction effects themselves directly. However, this result demonstrates that covariates that influence both PGSBMI 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.

Pearson correlations weighted by sample size between maximum R2 differences across strata, main effects of covariate on log(BMI), and PGS-covariate interaction effects on log(BMI). Main effect units are in standard deviations, interaction effect units are in PGS standard deviations multiplied by covariate standard deviations. Only continuous variables are plotted and modeled. GERA was excluded due to slightly different phenotype definitions.

Effects of machine learning approaches on predictive performance

Given evidence of PGS-covariate dependence, we aimed to assess increases in R2 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 PGSBMI), 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 R2 increased up to 67% (mean 23%) across cohorts and ancestries (Figure 5, Supplemental Table 7). The inclusion of interaction terms increased cross-validated R2 up to 12% (mean 3.9%) when comparing LASSO including interaction terms to LASSO with main effects only.

Model R2 across cohorts and ancestries using age and gender as covariates (along with PGSBMI and PCs 1-5). Across all cohorts and ancestries, LASSO with PGS-age and PGS-gender interaction terms had better average 10-fold cross-validation R2 than LASSO without interaction terms, while neural networks outperformed LASSO models.

We then modeled all available covariates for each cohort and did similar comparisons. Cross-validated R2 increased 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 suggest that machine learning methods can improve model R2 that include PGSBMI as variables beyond including interaction terms in linear models, even when variable selection is performed using LASSO.

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, 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 (PGSGxAge) increased test R2 to 0.0771 from 0.0715 (7.8% relative increase) compared to the best performing main effect PGS (Figure 6, Supplemental Table 9) – age-stratified PGS had much lower performance than both other strategies (unsurprising given reasons previously mentioned). Including a PGSGxAge-age interaction term only marginally increased R2 (.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 (Discussion), but we consider their investigation outside the scope of this study.

PGS R2 based on three sets of GWAS. “Main effects” were from a typical main effect GWAS, “GxAge” effects were from a GWAS with a SNP-age interaction term, and “Age stratified” GWAS had main effects only but were conducted in four age quartiles. PGS R2 was evaluated using two models: one with main effects only, and one with an additional PGS*Age interaction term.

Discussion

Across four large-scale cohorts of diverse ancestry, we uncover replicable effects of covariates on both performance and effects of PGSBMI. When stratifying by quintiles of different covariates, we find certain covariates that have significant and replicable evidence for differences in PGSBMI R2, with R2 being nearly double between top and bottom performing quintiles for covariates with the largest differences. When testing covariates for PGS-covariate interaction effects, we also find covariates with significant interaction effects, where, for the largest effect covariates, each standard deviation change affects PGSBMI effect by nearly 20%. Across analyses, we find 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. We also observed that interaction effects and R2 differences were strongly correlated, with main effects also being correlated with interaction effects and R2 differences, suggesting that covariates with the largest interaction effects also contribute to the largest R2 differences, with simple main effects also being predictive of expected differences in R2 and interaction effects. While in this study we stratified by quintile, similar trends of performance differences would be apparent in smaller percentiles. This observation is relevant to current research and clinical use of PGS, as individuals above a percentile cutoff are designated high risk (40). These analyses imply that individuals most at-risk for obesity should have both disproportionately higher predicted BMI and increased BMI prediction performance compared to the general population – more work is required to understand why, and if similar trends can be observed in other traits relevant to human health. We also employ machine learning methods for prediction of BMI with models that include PGSBMI and demonstrate that these methods outperform regularized linear regression models that include interaction effects. Finally, we employ 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.

Future work may include replicating these analyses across additional traits, and 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, the method of including a PGS for the covariate to better measure its environmental effect is potentially worth exploring further, which 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. Finally, incorporation of additional SNP-covariate interactions (e.g., SNP-gender) would also likely further improve prediction performance, although any SNP selection/adjustment procedures may further complicated by additional interaction terms.

We propose several ideas on why BMI-associated covariates have larger interaction effects and impact on R2 for PGSBMI. 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. Secondly, PGS may more readily explain more extreme phenotypes because they incorporate large effect variants associated with very high weight or height (41). For BMI, the distribution is skewed toward the right i.e., effects in trait-increasing variants may have a larger potential to affect trait variation compared to trait-reducing variants. This explanation would likely be better suited to traits with left-tailed distributions 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. In order to achieve accurate predictions as well as promote equity, further research is required regarding PGS methods and assessments. This research provides firm evidence on the context-specific nature of PGS and the incorporation of covariate interaction effects for predicting BMI and, apart from promoting equitable use of PGS across ancestries and cohorts, demonstrates ways to improve PGS performance incorporating context-specific effects.

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 (17) 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 < .250) were retained, and variants were filtered at INFO > .30 and minor allele frequency (MAF) > .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 (18) “-king-cutoff .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 (NEUR=376,729, NAFR=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 > .250) such that all were unrelated. 35,064 individuals (NEUR=31,961, NAFR=3,103) were used for analysis.

GERA

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

PMBB

Ancestry inference and relatedness inference have been described elsewhere (20). Individuals were split into European and African ancestry cohorts, and related individuals were removed at pi-hat > .250. 36,046 individuals (NEUR=26,372, NAFR=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,2123). 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 (PGSBMI) were constructed using PRS-CSx (24), using GWAS summary statistics for individuals of European (25), African (26), and East Asian (27) 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 8). The final PGSBMI per ancestry and cohort was calculated by regressing log(BMI) on the PGSBMI per ancestry without covariates – the combined, predicted value was used as a single PGSBMI 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 PGSBMI were used in the final models, as it had higher R2 than using the combined PGSBMI (data not shown).

PGSBMI 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 PGSBMI R2 was calculated by taking the difference in R2 between:

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

PGSBMI interaction modelling

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

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

Correlation between R2 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 (32)). We then estimated the correlation between main effects, interaction effects, and maximum R2 differences across all cohorts and ancestries weighting by sample size, analyzing quantitative and binary variables separately.

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). PGSBMI 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. ‘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 R2 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 [.01, .0001], type of learning rate [constant, inverse scaling], power t [.4, .6], momentum [.80, 1.0], batch size [32, 256], and number of hidden layers [1, 2].

GxAge PGSBMI 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 + PCs1-5

  2. BMI ∼ SNP + Age*SNP + Age + Sex + PCs1-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 (.50, .10, …, 5.0×10−5, 1.0×10−5) and remaining settings as default in Plink 1.9. For 2), GxAge PGSBMI 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:

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

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

For each of the three analyses, the parameters and models resulting in the best performing PGS (highest incremental R2, 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), R2 between all predicted values and true values was calculated together.

URLs

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

Data Availability statement

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.

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. 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).

Supplemental Figure

PGS-covariate interaction term p-values in UKBB EUR, with and without including the covariate PGS in the model – the mean -log10(p) is reduced from 18.0899 to 14.97072 with their inclusions. Note age and sex PGS were not calculated, and their interaction p-values are excluded from this figure.