Harnessing the power of polygenic risk scores to predict type 2 diabetes and its subtypes in a high-risk population of British Pakistanis and Bangladeshis in a routine healthcare setting

Study population

Genes & Health (G&H) recruits British Pakistani and British Bangladeshi (BPB) people aged 16 years and above, predominantly in community and primary care settings. We used the February 2020 data release(20) which comprises electronic health record (EHR) data and genotype data from the Illumina Infinium Global Screening Array V3 Chip on 22,944 participants(15). We applied specific inclusion and exclusion criteria to the G&H population for each analysis, described below. Descriptions of quality control and imputation of genotype data are provided in Supplementary File 1.

Ascertainment of conditions using clinical coding

Type 2 diabetes, gestational diabetes mellitus and associated complications were ascertained using clinical codes extracted from primary and secondary care EHR, presented in full in Table S1, developed using widely-used clinical coding resources (21). We excluded individuals with clinical codes of conditions causing secondary or monogenic diabetes.

Transfer of previously-identified GWAS loci to the study population

We obtained previously identified genome-wide association study (GWAS) loci associated with type 2 diabetes in people of European-ancestry (EUR) using Vujkovic et al. (22). We performed a GWAS of type 2 diabetes in G&H, and assessed if previously-identified loci were reproducible in G&H at P-value <0.05 (Supplementary File 1) using previously-established methods (6). The expected power for replication was estimated, assuming the same effect size as in the EUR discovery sample, and accounting for the allele frequency and sample size in G&H. Even though the same locus and target gene may affect type 2 diabetes risk in both populations, it is possible that this is caused by distinct causal variants in the region. Therefore, trans-ancestry colocalisation analysis was performed using TEColoc to assess whether a transferable locus shared the same causal variant between BPB and EUR populations, using UK Biobank (UKBB) EUR individuals for the latter (23). We applied the Popcorn algorithm to estimate the trans-ancestry genetic correlation between BPB and UKBB EUR populations (24).

Polygenic risk score construction

We applied previously constructed scores from the PGS Catalog(25), principally developed in EUR populations, to participants in G&H (see Table S2). These were compared to PRS optimised within G&H using the clumping and p-value thresholding (C+T) method implemented in PRSice2 v2.2.11 (26,27) based on the largest published type 2 diabetes multi-ancestry GWAS to date(22,25) (Figure S1). For the C+T score, LD estimated from EUR samples (N=503) from the 1000 Genomes project was used for clumping (r²=0.1). We calculated multiple scores using various P-value thresholds. We selected the PRS with the highest predictive accuracy in an independent set of samples that were not included in the QDiabetes analysis described below i.e. prevalent cases with onset before the assessment date, and controls who had insufficient data on the assessment date to be included. The area under the receiver operating characteristic curve (AUC) was estimated with the R package “pROC”. Predictive accuracy of PRS was quantified by incremental AUC, or the gain in AUC when adding PRS to the reference model which accounted for participant age, gender, and 10 genetic principal components (PCs). The C+T PRS showed higher predictive accuracy than scores from the PGS Catalog, and was thus used in downstream analysis (Figure S1). A scaled PRS following the normal distribution with mean and median of 0 and standard deviation of 1 was constructed across BPB participants after regressing out the first ten PCs, allowing direct comparison between ancestry groups(28).

Development of an integrated risk tool to predict incident type 2 diabetes

We assessed whether PRS can enhance 10-year risk prediction of type 2 diabetes compared to QDiabetes, a validated, EHR-based risk prediction tool commonly used in UK primary care to estimate an individual’s 10-year risk of developing type 2 diabetes(17). There are three QDiabetes models. Model A provides estimates based on age, ethnicity, family history, comorbidities, and prescribed medications. Model B (which has the highest clinical predictive value) uses the same variables as model A, plus fasting plasma glucose (FPG), while model C is composed of model A plus HbA1c. Analysis was performed in 13,648 individuals aged 25-84 years with no prior history of type 2 diabetes as of 1st May 2010, no recorded value of HbA1c > 48 mmol / mol or FPG > 7.0mmol/L (Table S1). These inclusion criteria aligned with those used for development and validation of the QDiabetes model (17). QDiabetes scores were calculated on the assessment date of 1st May 2010 for each participant (see Table 1 for participant characteristics) on the basis of available clinical data required to run the model. For each individual, the QDiabetes risk estimate was combined with the PRS to form an integrated estimate of risk of developing type 2 diabetes over the next 10 years (Integrated Risk Tool, IRT), as previously described (1). 10-year risk of disease was classified as high (>10%) or low (<10%) in line with (1). Relative model performance was assessed using the concordance index (C-Index) and net reclassification index (NRI; the proportion of participants correctly reclassified as high or low risk by the IRT compared to QDiabetes, minus the proportion of participants incorrectly reclassified) in all individuals and in age-by-sex subgroups (1,7). In a separate analysis restricted to women with a history of GDM, we compared characteristics (including type 2 diabetes PRS) of those who subsequently did (n = 127) and did not (n = 175) develop type 2 diabetes. Further details of incident disease prediction are provided in Supplementary File 1.

Identification of type 2 diabetes subgroups and their association with PRS

We applied data-driven clustering techniques to define subgroups of individuals with distinct characteristics in two separate models using latent class analysis (Stata V16.0) based on methodology used by (11–13). The optimal number of clusters was selected on the basis of information criteria elbow plots (Figure S2) and was found to be five for all model iterations. The first model (Clinical Model) was based on five clinical variables at the time of type 2 diabetes diagnosis: age, body mass index (BMI), HbA1c, serum triglycerides and serum high-density lipoprotein (HDL). Analysis was restricted to 4,266 participants for whom values for these variables were available within 1 year of type 2 diabetes diagnosis. Differences in PRS between clusters were compared using analysis of variance (ANOVA) (RStudio V1.1.413). In the second model (Integrated Model) we repeated clustering using age, BMI, HbA1c, serum triglycerides and PRS as covariates in 5,905 individuals with available data within 12 months of diagnosis. HDL was omitted from the Integrated Model as its inclusion did not alter cluster characteristics, but did reduce sample size. For both models, progression to micro- and macrovascular complications (defined in Table S1) was compared between clusters with Cox proportional hazard models controlling for hypertension, statin use, and serum total cholesterol using the R package “survival” (V2.3-11) in individuals free from micro- and macrovascular complications at the time of type 2 diabetes diagnosis.

Genetics of type 2 diabetes in BPB versus European individuals

We first investigated the extent to which genetic risk for type 2 diabetes was shared between G&H (BPB - British Pakistani and Bangladeshis) and EUR individuals. Results are described in more detail in Supplementary File 1. In brief, we found that trans-ethnic genetic correlation (i.e. the correlation of causal-variant effect sizes) between G&H and UKBB EUR individuals was significantly lower than 1 (r_g = 0.68, standard error = 0.15, p-value (for the null hypothesis that r_g = 1) = 0.03). Most of the 338 genetic loci ascertained in EUR populations which were well-powered to replicate in G&H did so (76 (22.5%) replicated versus 28.1% expected; binomial p=0.01; Table S3). Of 27 replicated loci assessed, 9 (33%) had significant evidence of shared causal variants (Table S4). For example, we observed shared causal variants at the TCF7L2 locus, one of the known loci with the strongest association with type 2 diabetes, and at the KCNJ11 locus, which is the target gene for drugs such as Glyburide (29,30).

Performance of an integrated risk tool to predict incident type 2 diabetes

We constructed a PRS for type 2 diabetes using the C+T method, based on multi-ancestry GWAS meta-analysis (22). This had an odds ratio (OR) per standard deviation (SD) of 1.57 (95% confidence interval [CI]: 1.50–1.65) and an incremental AUC of 0.032 (95% CI: 0.026, 0.039) in G&H. There was no correlation between PRS and QDiabetes scores (Pearson’s coefficients -0.03, 0.08, 0.13 for QDiabetes models A, B, and C respectively; associated p-values 0.31, 0.18, 0.16). Compared to QDiabetes alone, the integrated risk tool (IRT) combining QDiabetes model A (multiple clinical risk factors) with PRS improved the 10-year prediction of type 2 diabetes as assessed by net reclassification index (NRI); NRI 3.22% (95% CI: 2.00 - 4.38%) in 13,648 included individuals (Figure 1A). The IRT C-index, a goodness-of-fit metric approximating the area under the receiver-operator curve, was superior to QDiabetes model A for participants aged less than 40 years (p = 0.002) (Figure 2A / Table S5), but not for people aged over 40 years. Enhancement of type 2 diabetes prediction with the IRT persisted but was attenuated when QDiabetes model B (model A plus fasting plasma glucose) was used: NRI 0.80% (95% CI 0.21 - 1.42%). The IRT did not improve on the QDiabetes score with model C (model A plus HbA1c) (NRI 0.20%; 95% CI: -0.09 - 0.44%) (Figure 1B and 1C). There were no observed differences in C Index between QDiabetes and IRT in models B and C (Figures 2A and 2B / Table S5).

• Download figure
• Open in new tab

Figure 1:

Net Reclassification Index (NRI) and 95% confidence intervals comparing 10 year prediction of type 2 diabetes between QDiabetes models A (n = 13,648), B (n = 4,334), and C (n = 864) and integrated risk tools combining each QDiabetes model with type 2 diabetes polygenic risk score. NRI is presented overall and for age-by-sex subgroups. The NRI for individuals who subsequently develop type 2 diabetes (cases) and who do not develop type 2 diabetes (non-cases) are presented alongside the NRI for cases and non-cases combined (all individuals.

• Download figure
• Open in new tab

Figure 2:

Comparison of C Index calculated from Cox proportional hazard models for QDiabetes models A, B, and C and integrated risk tools combining each QDiabetes model with individual type 2 diabetes polygenic risk score. C index is presented for all participants, and for age-by-sex subgroups.

Across all IRT models (for QDiabetes model A, B and C), NRI was higher in participants aged less than 40 years who subsequently developed type 2 diabetes (cases) than those who did not (non-cases), implying improved ability to correctly classify younger individuals as high risk. The converse pattern was seen in individuals aged greater than 40 years, implying enhanced ability to correctly classify older individuals as low risk. The clinical features of individuals whose risk was reclassified are presented in Table S6. The number of individuals reclassified in each model are presented by age band and sex in Table S7.

Prediction of type 2 diabetes after gestational diabetes using PRS

In women who had a history of gestational diabetes (GDM), Type 2 diabetes PRS was higher in those who had subsequently developed type 2 diabetes, compared to those who had not (mean value (SD) 0.408 (0.93) vs 0.140 (0.98), p = 0.02) (Figure 3). Among individuals without type 2 diabetes (non-cases), women with a history of GDM displayed a higher mean PRS than age- and BMI-matched females without GDM, as well as males (p = 0.003 and p=0.001, respectively) (Figure 3). This could be because these individuals have increased propensity to develop type 2 diabetes in the future, but have not yet developed it, or because they have undiagnosed type 2 diabetes, for example, through missed screening tests post-delivery. In multivariate survival analysis restricted to women with a history of GDM, type 2 diabetes PRS was independently associated with development of type 2 diabetes after adjustment for (1) QDiabetes score model A (hazard ratio per SD of PRS (HR) 1.23, 95% CI 1.05 - 1.42; p = 0.028;) and (2) established risk factors for type 2 diabetes (HR 1.37, 95% CI 1.17 - 1.57, p = 0.002) (Table S8). Similar findings were observed in age- and BMI-matched male and female (without a history of GDM) controls (Table S8), indicating that the utility of the type 2 diabetes PRS for risk prediction does not appear to differ between women with versus without GDM.

• Download figure
• Open in new tab

Figure 3:

Mean type 2 diabetes polygenic risk score in women with gestational diabetes mellitus (GDM) who subsequently did and did not develop type 2 diabetes, compared to age- and BMI-matched control groups of females without GDM, and males. Data are presented as group mean with 95% confidence intervals according to sex-by-GDM subgroup, stratified by subsequent type 2 diabetes status.

Type 2 diabetes subgroup identification using clustering

In our Clinical Model, we identified type 2 diabetes subgroups by clustering age, BMI, HBA1c, HDL and triglycerides at the time of diagnosis (Figure 4A/Table 2). We replicated some previously identified clusters: Mild Obesity-related Diabetes (MOD), Mild Age-Related Diabetes (MARD), Severe Insulin Resistant Diabetes (SIRD) and Mild Diabetes (MD), and use their previous nomenclature (11–13,31). Our MD cluster may include individuals who are insulin deficient as it is characterised by low serum triglycerides (1.21 (0.43) mmol/L) in addition to high HDL (2.07 (0.25) mmol/L), but in the absence of clinical measures of insulin secretion such as C-peptide, this cannot be confirmed. Our cluster with the largest membership was undifferentiated by clinical features; we therefore call it Clinically-Undifferentiated Diabetes (CUD). There were no statistically significant differences in rates of progression to micro- or macrovascular complications between clusters in the Clinical model. However, we found that the PRS differed significantly between clusters(ANOVA p=0.002) and that CUD was characterised by the highest observed PRS (0.34 (0.97)). PRS was lowest in MOD and MARD (mean PRS (0.07 (1.02) and 0.23 (0.97)), respectively.

• Download figure
• Open in new tab

Figure 4:

Radar plots showing cluster mean valuess of the variables used to derive five clusters of individuals with type 2 diabetes at the time of diagnosis in two separate models. In panel A (Clinical Model), latent class analysis was performed using age, BMI (body mass index), HbA1c, serum triglycerides, and high density lipoprotein (HDL). In panel B (integrated model), clustering was repeated using type 2 diabetes polygenic risk score (PRS) in place of HDL. The centre and edge of each polygon represents the minimum and maximum mean cluster values for each variable, respectively. All polygons represent the same scale; all scales are linear. Mean values for each variable within each cluster are represented by coloured dots. For example, mean BMI is highest in cluster 1.

Given these results, we next explored whether inclusion of PRS as a covariate in the clustering model could further refine the subgroups (henceforth termed “Integrated model”). We included PRS in the latent class analysis along with age, BMI, HbA1c, and triglycerides (Figure 4B/Table 3). We identified broadly similar clusters, although with some additional differentiation not seen in the Clinical Model: MOD and MARD were maintained, but in the latter, mean age was notably higher than in the Clinical Model (61.6 (7.6) versus 50.9 (12.2) years). In this analysis, two clusters representing insulin resistant diabetes were identified, differentiated by both the severity of hypertriglyceridaemia and PRS, and therefore we call these SIRD and Insulin Resistant Diabetes (IRD). The MD cluster was not identified in this analysis. We again identified a prevalent clinically-undifferentiated cluster that was associated with a high PRS (mean PRS 0.48 (0.95)), and we therefore rename this cluster Clinically-Undifferentiated High Polygenic Susceptibility Diabetes (CUPS). In Cox proportional hazard models (Figure S3), macrovascular complication rates were higher in MARD and MOD (p<0.001 and p=0.009, respectively); microvascular rates were highest in IRD and MARD (p<0.001 and p=0.002).

Strengths and limitations

Strengths of this study include its assessment of an ethnic minority population with high levels of type 2 diabetes and associated morbidity and mortality that has been under-represented in research studies to date; its use of high quality longitudinal data from a relatively unselected sample and linkage to genetic data. Our genetic analysis is robust and we build a strong justification for the construction and evaluation of ancestry-specific PRSs. The main weaknesses of our work include our need to impute missing data not present in health records and the lack of a replication cohort to externally validate findings. Our analysis of the progression to type 2 diabetes after GDM was limited by the low uptake of postpartum diabetes screening that may have resulted in underdiagnosis of type 2 diabetes and ascertainment bias. Across all of our analyses, it is likely that some individuals coded as having type 2 diabetes actually have type 1 diabetes or rare monogenic forms of diabetes, although the absolute number of these and impact on overall findings is likely to be very small.

Conclusions

We highlight the need for ancestry-specific genetic tools to be developed for diverse populations to elucidate aetiological differences. Our PRS, optimised in British South Asians, enhances prediction of incident type 2 diabetes when combined with an established clinical risk tool, compared to using the clinical tool alone, and particularly in young people. Additionally, the PRS has value in predicting the onset of type 2 diabetes in a specific high risk group, women affected by gestational diabetes. These findings could aid the personalisation of care of people at risk of type 2 diabetes. The PRS also adds a unique dimension to the identification of diabetes subgroups at diagnosis, and their association with future complications, that may be valuable for stratification of care in the future. Our work provides important insight into the genetic risk in an ethnic group underrepresented in research but disproportionately affected by type 2 diabetes, and has significant potential to be translated into clinical practice.