Multi-ethnic polygenic risk scores improve risk prediction in diverse populations

TL;DR

A new method called MultiPRS improves disease risk prediction in non-European populations by combining data from European and target population samples.

Contribution

MultiPRS introduces a novel approach that combines large European and small target population data for better risk prediction.

Findings

MultiPRS achieved over 70% relative improvement in predicting type 2 diabetes in a Latino cohort.

Improvements were consistent with simulation results and depend on ancestry-adjusted coefficients.

The method reduces the accuracy gap in genetic risk prediction between European and non-European populations.

Abstract

Methods for genetic risk prediction have been widely investigated in recent years. However, most available training data involves European samples, and it is currently unclear how to accurately predict disease risk in other populations. Previous studies have used either training data from European samples in large sample size or training data from the target population in small sample size, but not both. Here, we introduce a multi-ethnic polygenic risk score that combines training data from European samples and training data from the target population. We applied this approach to predict type 2 diabetes (T2D) in a Latino cohort using both publicly available European summary statistics in large sample size and Latino training data in small sample size. We attained a >70% relative improvement in prediction accuracy (from R2=0.027 to R2=0.047) compared to methods that use only one source of training data, consistent with large relative improvements in simulations. We observed a systematically lower load of T2D risk alleles in Latino individuals with more European ancestry, which could be explained by polygenic selection in ancestral European and/or Native American populations. Application of our approach to predict T2D in a South Asian UK Biobank cohort attained a >70% relative improvement in prediction accuracy, and application to predict height in an African UK Biobank cohort attained a 30% relative improvement. Our work reduces the gap in polygenic risk prediction accuracy between European and non-European target populations.

The use of genetic information to predict disease risk is of great interest because of its potential clinical application. Prediction is performed via the construction of polygenic risk scores, which separate individuals into different risk categories. Polygenic risk scores can also be applied to improve our understanding of the genetic architecture of complex diseases. The ideal training data set would be a large cohort from the same population as the target sample, but this is generally unavailable for non-European populations. Thus, we propose a summary statistics based polygenic risk score that leverages both a large European training sample and a training sample from the same population as the target population. This approach produces a substantial relative improvement in prediction accuracy compared to methods that use a single training population when applied to predict type 2 diabetes in a Latino cohort, consistent with simulation results. We observed similar relative improvements in applications to predict type 2 diabetes in a South Asian cohort and height in an African cohort.

Figures

Accuracy of 5 prediction methods in simulations using subsets of chromosomes.We report prediction accuracies for each of the 5 main prediction methods as a function of M/Msim, where M=232,629 is the total number of SNPs and Msim is the actual number of SNPS used in each simulation: 232,629 (all chromosomes), 68,188 (chromosomes 1-4), 38,412 (chromosomes 1-2), and 19,087 (chromosome 1). Numerical results are provided in S5 Table.

Accuracy of 5 prediction methods in analyses of type 2 diabetes in a Latino cohort as a function of P-value thresholds.We report prediction accuracies for each of the 5 main prediction methods as a function of (a) EUR P-value threshold, where applicable (with optimized LAT P-value threshold, where applicable) and (b) LAT P-value threshold, where applicable (with optimized EUR P-value threshold, where applicable). Numerical results are provided in S13 Table and S14 Table.

Tables

Full text

Introduction

Genetic risk prediction is an important and widely investigated topic because of its potential clinical application as well as its application to better understand the genetic architecture of complex traits [1]. Many polygenic risk prediction methods have been developed and applied to complex traits. These include polygenic risk scores (PRS) [2–9], which use summary association statistics as training data, and Best Linear Unbiased Predictor (BLUP) methods and their extensions [10–17], which require individual-level genotype and phenotype data.

However, all of these methods are inadequate for polygenic risk prediction in non-European populations, because they consider training data from only a single population. Existing training data sets have much larger sample sizes in European populations, but the use of European training data for polygenic risk prediction in non-European populations reduces prediction accuracy, due to different patterns of linkage disequilibrium (LD) in different populations [2,8,18,19]. For example, ref. [8] reported a relative decrease of 53-89% in schizophrenia risk prediction accuracy in Japanese and African-American populations compared to Europeans when applying PRS methods using European training data. An alternative is to use training data from the same population as the target population, but this would generally imply a much lower sample size, reducing prediction accuracy.

To tackle this problem, we developed an approach that combines PRS based on European training data with PRS based on training data from the target population. The method takes advantage of both the accuracy that can be achieved with large training samples [4,5] and the accuracy that can be achieved with training data containing the same LD patterns as the target population. In simulations and application to predict type 2 diabetes (T2D) in Latino target samples in the SIGMA T2D data set [20], we attained a >70% relative improvement in prediction accuracy (from R^2^=0.027 to R^2^=0.047) compared to methods that use only one source of training data. We also obtained a >70% relative improvement in an analysis to predict T2D in a South Asian UK Biobank cohort, and a 30% relative improvement in an analysis to predict height in an African UK Biobank cohort.

Materials and Methods

Polygenic risk score using a single training population

Polygenic risk scores are constructed using SNP effect sizes estimated from genome-wide association studies, which perform marginal regression of the phenotype of interest on each SNP in turn. Explicitly, for continuous traits, we estimate effect sizes [eqn] (where i = 1,…,M indexes genetic markers) using the model y = b0 + b_i_g_i_ + b_PC_PC + ε, where g_i_ denotes genotypes at marker i, PC denotes one or more principal components used to adjust for ancestry, and ε denotes environmental noise. For binary traits, we use the analogous logistic model logit[P(y=1)] = b0 + b_i_g_i_ + b_PC_PC + ε.

Given a vector of estimated effect sizes [eqn] from a genome-wide association study performed on a set of training samples, the polygenic risk score [2] (PRS) for a target individual with genotypes g_i_ is defined as [eqn] In practice, rather than computing the PRS using estimated effect sizes for all available genetic markers, the PRS is computed on a subset of genetic markers obtained via informed LD-pruning [3] (also known as LD-clumping) followed by P-value thresholding [2]. Specifically, this “pruning + thresholding” strategy has two parameters, R_LD_^2^ and P_T_, and proceeds as follows. First, we prune the SNPs based on a pairwise threshold R_LD_^2^, removing the less significant SNP in each pair (using PLINK; see Web Resources). Second, we restrict to SNPs with an association P-value below the significance threshold P_T_.

The parameters R_LD_^2^ and P_T_ are commonly tuned using on validation data to optimize prediction accuracy [2,3]. While in theory this procedure is susceptible to overfitting, in practice, validation sample sizes are typically large, and R_LD_^2^ and P_T_ are selected from a small discrete set of parameter choices, so overfitting is considered to have a negligible effect. Accordingly, in this work, we consider R_LD_^2^ ∈ {0.1, 0.2, 0.5, 0.8} and P_T_ ∈ {1.0, 0.8, 0.5, 0.4, 0.3, 0.2, 0.1, 0.08, 0.05, 0.02, 0.01, 10^−3^, 10^−4^, 10^−5^, 10^−6^, 10^−7^, 10^−8^}, and we always report results corresponding to the best choices of these parameters. In all of our primary analyses involving two training populations (see below), values of R_LD_^2^ and P_T_ were optimized based only on PRS in a single training population, to ensure that PRS using two training populations did not gain any relative advantage from the optimization of these parameters.

In this work, we specifically consider PRS built using European (EUR), Latino (LAT), South Asian (SAS), or African (AFR) training samples. We use the notation PRS_EUR_ to denote PRS built using European samples, and analogously for the other populations.

Polygenic risk score using two training populations

Given a pair of polygenic risk scores computed as above using two distinct training populations, we define the multi-ethnic PRS with mixing weights α_1_ and α_2_ as the linear combination of the two PRS with these weights: e.g., for EUR and LAT, we define PRS_EUR+LAT_ = α_1_PRS_EUR_ + α_2_PRS_LAT_. We employ two different approaches to avoid overfitting. In our primary analyses, we estimate mixing weights α_1_ and α_2_ using validation data and compute adjusted R^2^ to account for the additional degree of freedom. In our secondary analyses, we estimate mixing weights α_1_ and α_2_ using cross-validation (see Assessment of methods below).

Polygenic risk score using one or two training populations and genetic ancestry

We further define polygenic risk scores that include an ancestry predictor, namely, the top principal component in a given data set. (We considered only the top PC in each data set that we analyzed, because lower PCs had a squared correlation with phenotype lower than 0.005 in each case.) We define a polygenic risk score LAT+ANC with mixing weights α_1_ and α_2_ as PRS_LAT+ANC_ = α_1_PRS_LAT_ + α_2_ PC, and we define a polygenic risk score EUR+LAT+ANC with mixing weights α_1_, α_2_ and α_3_ as PRS_EUR+LAT+ANC_ = α_1_PRS_EUR_ + α_2_PRS_LAT_ + α_3_PC. As above, we employ two different approaches to avoid overfitting: in our primary analyses, we estimate mixing weights using validation data and compute adjusted R^2^; in our secondary analyses, we estimate mixing weights using cross-validation.

Assessment of methods

We assessed the accuracy of polygenic risk scores in validation samples (independent from samples used to estimate effect sizes). We used adjusted R^2^ as the accuracy metric for continuous traits and liability-scale adjusted R^2^ (ref. [21]) for binary traits. Adjusted R^2^ is defined as [eqn] is the number of PRS or ANC components in the mixture, n is the number of validation samples, and [eqn] is the raw (unadjusted) R^2^. The adjusted R^2^ metric roughly corrects for increased model complexity in multi-component PRS, so in our primary analyses, we report accuracy as adjusted R^2^ using best-fit mixing weights [eqn] estimated using the validation data.

To verify that this metric provides robust model comparisons, we also performed auxiliary analyses in which we used 10-fold cross-validation: specifically, for each left-out fold in turn, we estimated mixing weights using the other 9 folds and evaluated adjusted R^2^ for PRS computed using these weights on the left-out fold. We then computed average adjusted R^2^ across the 10 folds. (When analyzing data from an unbalanced case-control study with #cases << #controls, we used stratified 10-fold cross-validation, selecting the folds such that each fold had the same case-control ratio; this applies only to the South Asian UK Biobank T2D analysis.)

Finally, for analyses in which we needed to use samples from the same cohort for both building PRS (i.e., estimating effect sizes [eqn]) and validation, we also used cross-validation. In our primary analyses, we employed 10-fold cross-validation, using 90% of the cohort to estimate [eqn] and the remaining 10% of the cohort to validate predictions (using the adjusted R^2^ metric with best-fit mixture weights [eqn]). In our secondary analyses, we employed 10×9-fold cross-validation, in which 90% of the cohort was used to estimate both [eqn] and [eqn] and the remaining 10% of the cohort was used to validate predictions. To estimate [eqn], we iteratively split the 90% set of training samples into an 80% training-training set and a 10% training-test set; we estimated [eqn] in the 80% training-training set and computed a PRS for the 10% training-test set for each of the 9 training-test folds, and we then performed a single regression of phenotype against each PRS across the entire 90% set of training samples to estimate [eqn]. Finally, we reestimated [eqn] for the final test prediction using the entire 90% set of training samples.

Simulations

We simulated quantitative phenotypes using real genotypes from European (WTCCC2) and Latino (SIGMA) data sets (see below). We fixed the proportion of causal markers at 1% and fixed SNP-heritability h_g_^2^ at 0.5, and sampled normalized effect sizes β_i_ from a normal distribution with variance equal to h_g_^2^ divided by the number of causal markers. We calculated per-allele effect sizes b_i_ as [eqn] where p_i_ is the minor allele frequency of SNP i in the European data set. We simulated phenotypes as [eqn]

In our primary simulations, we discarded the causal SNPs and used only the non-causal SNPs as input to the prediction methods (i.e. we simulated untyped causal SNPs, which we believe to be realistic). As an alternative, we also considered simulations in which we included the causal SNPs as input to the prediction methods (i.e., a scenario in which causal SNPs are typed).

We also performed simulations in which Latino phenotypes were explicitly correlated to ancestry (population stratification). In these simulations, we added a constant multiple of PC1 (representing European vs. Native American ancestry, with positive values representing higher European ancestry) to the Latino phenotypes such that the correlation between phenotype and PC1 was equal to -0.11, which is the correlation between the T2D phenotype and PC1 in the SIGMA data set.

We performed simulations under 4 different scenarios: (i) using all chromosomes, (ii) using chromosomes 1-4, (iii) using chromosomes 1-2, and (iv) using chromosome 1 only. The motivation for performing simulations with a subset of chromosomes was to increase N/M, extrapolating to performance at larger sample sizes, as in previous work [8].

Simulation data sets: WTCCC2 and SIGMA

Our simulations used real genotypes from the WTCCC2 and SIGMA data sets. The WTCCC2 data set consists of 15,622 unrelated European samples from a multiple sclerosis study genotyped at 360,557 SNPs after QC [22,23] (see Web Resources). The SIGMA data set consists of 8,214 unrelated Latino samples genotyped at 2,440,134 SNPs after QC [20] (see Web Resources). We restricted our simulations to 232,629 SNPs present in both data sets (with matched reference and variant alleles) after removing A/T and C/G SNPs to eliminate potential strand ambiguity.

Latino type 2 diabetes data sets: DIAGRAM, SIGMA and UK Biobank

Our analyses of type 2 diabetes in Latinos used summary association statistics from the DIAGRAM data set and genotypes and phenotypes from the SIGMA data set. The DIAGRAM data set consists of 12,171 cases and 56,862 controls of European ancestry for which summary association statistics at 2,473,441 imputed SNPs are publicly available (see Web Resources) [24]. As noted above, the SIGMA data set consists of 8,214 unrelated Latino samples (3,848 type 2 diabetes cases and 4,366 controls) genotyped at 2,440,134 SNPs after QC. SIGMA association statistics were computed with adjustment for 2 PCs, as in ref. [20]. We restricted our analyses of type 2 diabetes to 776,374 SNPs present in both data sets (with matched reference and variant alleles) after removing A/T and C/G SNPs to eliminate potential strand ambiguity. For the SIGMA data set, we used the top 2 PCs as computed in ref. [20]. We also performed an analysis of type 2 diabetes using imputed genotypes from the SIGMA T2D data set [20], restricting to 2,062,617 SNPs present in both data sets (with matched reference and variant alleles) after removing A/T and C/G SNPs to eliminate potential strand ambiguity.

We performed a secondary analysis using 113,851 British samples from UK Biobank [25] (see Web Resources) as European training data (5,198 type 2 diabetes cases and 108,653 controls). UK Biobank association statistics were computed with adjustment for 10 PCs [25], estimated using FastPCA [26] (see Web Resources). We computed summary statistics for 608,878 genotyped SNPs from UK Biobank after removing A/T and C/G SNPs to eliminate potential strand ambiguity. We analyzed 187,142 SNPs present in the SIGMA and UK Biobank data sets. We defined type 2 diabetes cases in UK Biobank as “any diabetes” with “age of diagnosis > 30”. We note that the p-values at two top type 1 diabetes (T1D) loci (rs2476601, rs9268645) were only nominally significant (p~0.05) for this T2D phenotype, indicating low contamination with T1D cases.

South Asian type 2 diabetes data sets: DIAGRAM, SAT2D and UK Biobank

Our analysis of type 2 diabetes in South Asians used European summary association statistics from the DIAGRAM data set (described above), South Asian summary statistics data from the South Asian Type 2 Diabetes (SAT2D) Consortium [27], and South Asian genotypes and phenotypes from UK Biobank (see Web Resources) as test data. The SAT2D data set consists of 5,561 South Asian type 2 diabetes cases and 14,458 South Asian controls for which we summary statistics for 2,646,472 imputed SNPs were available. The UK Biobank test data consists of 1,756 unrelated samples of South Asian ancestry (272 type 2 diabetes cases and 1,484 controls), genotyped at 608,878 SNPs after QC, with the following self-reported ethnicity distribution: 52 Bangladeshi, 1,301 Indian and 403 Pakistani. We analyzed 208,400 SNPs present in the DIAGRAM, SAT2D and UK Biobank data sets after removing A/T and C/G SNPs to eliminate potential strand ambiguity.

African height data sets: UK Biobank and N’Diaye et al

Our analyses of height in Africans used European summary association statistics from UK Biobank (see Web Resources), African summary statistics from N’Diaye et al. [28] and African genotypes and phenotypes from UK Biobank. European summary statistics from UK Biobank were computed using 113,660 British samples for which height phenotypes were available with adjustment for 10 PCs [25], estimated using FastPCA [26] (see Web Resources). The N’Diaye et al. [28] data set consists of 20,427 samples of African ancestry with summary association statistics at 3,254,125 imputed SNPs. The UK Biobank data set consists of 1,745 unrelated samples of African ancestry, genotyped at 608,878 SNPs after QC, with the following self-reported ethnicity distribution: 743 African, 1,002 Caribbean. We restricted our analysis to 232,182 SNPs present in the UK Biobank and N’Diaye et al. data sets after removing A/T and C/G SNPs to eliminate potential strand ambiguity.

Results

Simulations

We performed simulations using real genotypes and simulated phenotypes. We simulated continuous phenotypes under a non-infinitesimal model with 1% of markers chosen to be causal with the same effect size in all samples and SNP-heritability h_g_^2^ = 0.5 (see Methods); we report the average adjusted R^2^ and standard errors over 100 simulations. We used WTCCC2 [22,23] data (15,622 samples after QC; see Methods) as the European training data, and the SIGMA data [20] (8,214 samples) as the Latino training and validation data (with 10-fold cross-validation). We simulated phenotypes using the 232,629 SNPs present in both data sets and built predictions from these SNPs excluding the causal SNPs, modeling the causal SNPs as untyped (see Methods).

Prediction accuracies (adjusted R^2^) and optimal weights for the 5 main methods (EUR, LAT, LAT+ANC, EUR+LAT, EUR+LAT+ANC) are reported in Table 1. In each case, the best prediction accuracy was attained using LD-pruning threshold RLD^2^=0.8 (results using different LD-pruning thresholds are reported in S1 Table); the median value of the optimal P-value threshold P_T_ was equal to 0.01 for EUR and 0.05 for LAT. On average, the EUR method performed only 23% better than the LAT method, despite having twice as much training data. This reflects a tradeoff between the larger training sample size for EUR and the target-matched LD patterns for LAT. EUR+LAT attained 64%-101% relative improvements vs. EUR and LAT respectively (and used a slightly larger weight for EUR than for LAT), highlighting the advantages of incorporating multiple sources of training data. When including an ancestry predictor, EUR+LAT+ANC attained a 10% relative improvement vs. EUR+LAT (≥80% relative improvement vs. EUR or LAT), reflecting small genetic effects of ancestry on phenotype that can arise from random genetic drift between populations at causal markers (which is better-captured by ancestry components than by SNPs used in a PRS). Predictions using Latino effect sizes that were not adjusted for genetic ancestry (LAT_unadj_, EUR+LAT_unadj_, EUR+LAT_unadj_+ANC, as compared to LAT, EUR+LAT, EUR+LAT+ANC) were much less accurate (S2 Table), as in previous work [29]; this is consistent with the fact that LAT_unadj_ predictions were dominated by genetic ancestry (adjusted R^2^ = 0.37; S3 Table). We also observed a modest correlation (adjusted R^2^ = 0.025) between the EUR prediction and genetic ancestry (S3 Table), again reflecting small genetic effects of ancestry on phenotype that can arise from random genetic drift between populations at causal markers. The relative performance of the different prediction methods was similar in simulations in which phenotypes explicitly contained an ancestry term, representing environmentally-driven stratification (S4 Table).

Table 1.: Accuracy of 5 prediction methods in simulations.We report average adjusted R2 over 100 simulations for each of the 5 main prediction methods. We also report normalized weights, defined as the mixing weight (see Methods) multiplied by the standard deviation of the PRS.

We extrapolated the results in Table 1 to larger sample sizes by limiting the simulations to subsets of chromosomes, as in previous work [8] (Fig 1 and S5 Table). EUR+LAT+ANC was the best performing method in each of these experiments. We also performed simulations using predictions constructed using all SNPs including the causal SNPs (S1 Fig and S6 Table). In these experiments, EUR+LAT+ANC was once again the best performing method, and EUR performed much better than LAT, consistent with the larger training sample size for EUR and the fact that differential tagging of causal SNPs is of reduced importance when causal SNPs are typed.

Accuracy of 5 prediction methods in simulations using subsets of chromosomes.We report prediction accuracies for each of the 5 main prediction methods as a function of M/Msim, where M=232,629 is the total number of SNPs and Msim is the actual number of SNPS used in each simulation: 232,629 (all chromosomes), 68,188 (chromosomes 1-4), 38,412 (chromosomes 1-2), and 19,087 (chromosome 1). Numerical results are provided in S5 Table.

Analyses of type 2 diabetes in Latinos

We applied the same methods to predict T2D in Latino target samples from the SIGMA T2D data set. We used publicly available European summary statistics from DIAGRAM[24] (12,171 cases and 56,862 controls; effective sample size = 4/(1/Ncase + 1/Ncontrol) = 40,101) as European training data and SIGMA T2D genotypes and phenotypes [20] (3,848 cases and 4,366 controls; effective sample size = 8,181) as Latino training and validation data, employing 10-fold cross-validation.

Prediction accuracies (adjusted R^2^ on the liability scale [21], assuming 8% prevalence [3]) and optimal weights for the 5 main methods (EUR, LAT, LAT+ANC, EUR+LAT, EUR+LAT+ANC) are reported in Table 2 (other prediction metrics are reported in S7 Table). In each case, the best prediction accuracy was obtained using LD-pruning threshold RLD^2^=0.8 (results using different LD-pruning thresholds are reported in S8 Table); the value of the optimal P-value threshold PT was equal to 0.05 for EUR and 0.2 for LAT. EUR performed only 33% better than LAT despite the much larger training sample size, again reflecting a tradeoff between sample size and target-matched LD patterns. EUR+LAT attained 75%-133% relative improvements vs. EUR and LAT respectively (and used a slightly larger weight for EUR than for LAT), again highlighting the advantages of incorporating multiple sources of training data. Although adding an ancestry predictor to LAT produced a substantial improvement (LAT+ANC vs. LAT), adding an ancestry predictor to EUR+LAT produced an insignificant change in accuracy for EUR+LAT+ANC compared to EUR+LAT; this can be explained by the large negative correlation between the European PRS (EUR) and the proportion of European ancestry within Latino samples (R = −0.75; S9 Table), such that any predictor that includes EUR already includes effects of genetic ancestry. This correlation is far larger than analogous correlations due to random genetic drift in our simulations (S3 Table), suggesting that this systematically lower load of T2D risk alleles in Latino individuals with more European ancestry could be due to polygenic selection [30,31] in ancestral European and/or Native American populations; previous studies using top GWAS-associated SNPs have also reported continental differences in genetic risk for T2D [32,33]. We observed a similar correlation (R= −0.77) when using British UK Biobank type 2 diabetes samples as European training data (see Methods), confirming that this negative correlation is not caused by population stratification in DIAGRAM. As in our simulations, predictions using Latino effect sizes that were not adjusted for genetic ancestry (LAT_unadj_, EUR+LAT_unadj_, EUR+LAT_unadj_+ANC, as compared to LAT, EUR+LAT, EUR+LAT+ANC) were much less accurate (S10 Table), consistent with the fact that these predictions were dominated by genetic ancestry (S9 Table). We also computed predictions for each method using imputed SNPs from the SIGMA T2D data set; this did not improve prediction accuracy, but predicting using two training populations still achieved the highest accuracy (Supplementary Table 11).

Table 2.: Accuracy of 5 prediction methods in analyses of type 2 diabetes in a Latino cohort.We report adjusted R2 on the liability scale for each of the 5 main prediction methods. We obtained similar relative results using Nagelkerke R2, R2 on the observed scale and AUC (S7 Table). P-values are from likelihood ratio tests comparing models EUR and LAT to the null model, model LAT+ANC to LAT, model EUR+LAT to EUR, and EUR+LAT+ANC to EUR+LAT. For the EUR model we used RLD2=0.8 and PT=0.05 and for LAT we used RLD2=0.8 and PT=0.2. We also report normalized weights, defined as the mixing weight (see Methods) multiplied by the standard deviation of the PRS.

We investigated how the prediction accuracy of each method varied as a function of P-value thresholds, by varying either the EUR P-value threshold (Fig 2a and S12 Table) or the LAT P-value threshold (Fig 2b and S13 Table) between 10^-8^ and 1. In both cases, permissive P-value thresholds performed best, reflecting the relatively small sample sizes analyzed. However, the prediction accuracy of EUR+LAT+ANC was relatively stable, with prediction adjusted R^2^ > 0.037 across all EUR P-value thresholds (Fig 2a) and adjusted R^2^ > 0.033 across all LAT P-value thresholds (Fig 2b). In Fig 2a, we observe that as the EUR P-value threshold becomes more stringent, the difference in prediction accuracy between EUR+LAT+ANC and EUR+LAT increases, because EUR is less able to capture polygenic ancestry effects (see above).

Accuracy of 5 prediction methods in analyses of type 2 diabetes in a Latino cohort as a function of P-value thresholds.We report prediction accuracies for each of the 5 main prediction methods as a function of (a) EUR P-value threshold, where applicable (with optimized LAT P-value threshold, where applicable) and (b) LAT P-value threshold, where applicable (with optimized EUR P-value threshold, where applicable). Numerical results are provided in S13 Table and S14 Table.

In the above results (Table 2 and Fig 2), we allowed each prediction method to optimize its mixing weights via an in-sample fit in the target sample. This procedure could in principle be susceptible to overfitting [34,35]. We did not expect overfitting to be a concern given the small number of mixing weights optimized (at most 3) relative to the target sample size (8,181) and given our use of adjusted R^2^ as the evaluation metric, but to verify this expectation, we repeated our analyses using 10x9-fold cross-validation (see Methods). Methods that use two training populations remained much more accurate than single ancestry methods, as prediction accuracy decreased only very slightly (2-4% relative decrease vs. Table 2) for each method (S14 Table). These slight decreases are expected, since mixing weights optimized within 10x9 cross-validation are slightly suboptimal (due to reduced training data) and prediction accuracy is mildly sensitive to the choice of mixing weights (S2 Fig).

Analyses of type 2 diabetes in South Asians

We applied the same methods to predict T2D in South Asian target samples from the UK Biobank. We used publicly available European summary statistics from DIAGRAM (12,171 cases and 56,862 controls; effective sample size = 40,101) as European training data, South Asian summary statistics from SAT2D [27] (5,561 cases and 14,458 controls; effective sample size = 16,065) as South Asian training data, and UK Biobank genotypes and phenotypes (272 cases and 1,484 controls; effective sample size = 919) as South Asian validation data (see Methods).

Prediction accuracies (adjusted R^2^ on the liability scale [21], assuming sample prevalence 15%) and optimal weights for the 5 main methods (EUR, SAS, SAS+ANC, SAS+LAT, EUR+SAS+ANC) are reported in Table 3 (other prediction metrics are reported in S15 Table). In each case, the best prediction accuracy was obtained using LD-pruning threshold RLD^2^=0.8 (results using different LD-pruning thresholds are reported in S16 Table); the value of the optimal P-value threshold P_T_ was equal to 10^−3^ for EUR and 0.8 for SAS. EUR performed only 14% better than SAS despite the larger training sample size, again reflecting a tradeoff between sample size and target-matched LD patterns. EUR+SAS attained 72%-95% relative improvements vs. EUR and SAS respectively (and used a slightly larger weight for EUR than for SAS). Adding an ancestry predictor to EUR+SAS produced an insignificant change in accuracy for EUR+ SAS +ANC compared to EUR+SAS; we note a modest correlation between each prediction method and the proportion of European-related ancestry [36] within South Asian samples (see S17 Table). We repeated our analyses using stratified 10-fold cross-validation to estimate mixing weights (see Methods). We observed that methods that use two training populations continued to substantially outperform PRS using a single training population despite a decrease in prediction adjusted R^2^ (vs. Table 3) for each method, consistent with the limited sample size for estimating mixing weights (S18 Table).

Table 3.: Accuracy of 5 prediction methods in analyses of type 2 diabetes in a South Asian cohort.We report adjusted R2 on the liability scale for each of the 5 main prediction methods. We obtained similar relative results using Nagelkerke R2, R2 on the observed scale and AUC (S16 Table). P-values are from likelihood ratio tests comparing models EUR and SAS to the null model, model SAS+ANC to SAS, model EUR+SAS to EUR, and EUR+LAT+ANC to EUR+SAS. For the EUR model we used RLD2=0.8 and PT=10−3 and for SAS we used RLD2=0.8 and PT=0.8. We also report normalized weights, defined as the mixing weight (see Methods) multiplied by the standard deviation of the PRS.

Analyses of height in Africans

We applied the same methods to predict height in African target samples from the UK Biobank. We used European summary statistics from UK Biobank (113,660 samples; British ancestry only) as European training data, African summary statistics from ref. [28] (20,427 samples) as African training data, and African UK Biobank genotypes and phenotypes (1,745 samples) as African validation data.

Prediction accuracies (adjusted R^2^) and optimal weights for the 5 main methods (EUR, AFR, AFR+ANC, EUR+AFR, EUR+AFR+ANC) are reported in Table 4. For EUR and AFR, the best prediction accuracy was obtained using R_LD_^2^=0.2 and R_LD_^2^=0.8 respectively, thus we used these respective values of R_LD_^2^for EUR and AFR in each PRS in all primary analyses (results using different LD thresholds are reported in S19 Table); the value of the optimal P-value threshold P_T_ was equal to 10^−3^ for EUR and 0.05 for AFR. EUR performed much better than AFR, consistent with the far larger training sample size.

Table 4.: Accuracy of 5 prediction methods in analyses of height in an African cohort.We report adjusted R2 on the observed scale for each of the 5 main prediction methods. P-values are from likelihood ratio tests comparing models EUR and AFR to the null model, model AFR+ANC to AFR, model EUR+AFR to EUR, and EUR+LAT+ANC to EUR+AFR. For the EUR model we used RLD2=0.2 and PT=10−3 and for AFR we used RLD2=0.8 and PT=0.05. We also report normalized weights, defined as the mixing weight (see Methods) multiplied by the standard deviation of the PRS.

Nevertheless, EUR+AFR attained a 30% improvement vs. EUR (using a larger weight for EUR than for AFR). Adding an ancestry predictor to EUR+AFR produced an insignificant change in accuracy for EUR+AFR+ANC compared to EUR+AFR; we note a modest correlation between each prediction method and the proportion of European-related ancestry [36] within African samples (see S20 Table). We repeated our analyses using stratified 10-fold cross-validation to estimate mixing weights (see Methods). We observed that methods that use two training populations continued to substantially outperform PRS using a single training population despite a decrease in prediction adjusted R^2^ (vs. Table 4) for each method, consistent with the limited sample size for estimating mixing weights (S21 Table).

Discussion

We have shown that combining training data from European samples and training data from the target population attains a >70% relative improvement in prediction accuracy for type 2 diabetes in both Latino and South Asian cohorts compared to prediction methods that use training data from a single population. In addition, this approach attains 30% relative improvement in prediction accuracy for height in an African cohort. These relative improvements are robust to overfitting, consistent with simulations and reduce the documented gap in risk prediction accuracy between European and non-European target populations [2,8,18,19,37,38]. Intuitively, our approach leverages both large training sample sizes and training data with target-matched LD patterns. We note that the effects of differential tagging (or different causal effect sizes) in different populations can potentially be quantified using cross-population genetic correlation [39–41], and that leveraging data from a different population to improve predictions is a natural analogue to leveraging data from a correlated trait [14].

Despite these advantages, our work is subject to limitations and leaves several questions open for future exploration. First, although we have demonstrated large relative improvements in prediction accuracy, absolute prediction accuracies are currently not large enough to achieve clinical utility, which will require larger sample sizes [4,5]; our simulations suggest that multi-ethnic polygenic risk scores will continue to produce improvements at larger sample sizes (Fig 1). Second, while our focus here was on prediction without using individual-level training data, when such data is available it may be possible to attain higher prediction accuracy using methods that fit all markers simultaneously, such as Best Linear Unbiased Predictor (BLUP) methods and their extensions [10–17]. Third, our LDpred risk prediction method [8], which analyzes summary statistics in conjunction with LD information from a reference panel, is more accurate in European populations than the informed LD-pruning + P-value thresholding approach employed here; we did not employ LDpred due to the complexities of admixture-LD in analyses of admixed populations that explicitly model LD [42], but extending LDpred to handle these complexities could further improve accuracy. Fourth, we note that in our application to real phenotypes adding an ancestry predictor produced insignificant changes in prediction accuracy, primarily because ancestry effects are captured by the polygenic risk scores; adding an ancestry predictor only improves prediction when we use a stringent P-value threshold to build the polygenic risk score (Fig 2). Fifth, we have not considered here how to improve prediction accuracy in data sets with related individuals [15]. Finally, we focused our analyses on common variants, but future work may wish to consider rare variants as well.