Nuclear magnetic resonance-based metabolomics with machine learning for predicting progression from prediabetes to diabetes

Study design and participants

The UK Biobank is a large population-based prospective cohort study enrolling more than 500,000 community-dwelling adults from 22 assessment centers across the UK between 2006 and 2010^20,21. Participants completed touchscreen questionnaires and physical measurements and provided blood samples at baseline. The study was approved by the Northwest Multicenter Research Ethics Committee (REC reference for UK Biobank 11/NW/0382), and all participants provided informed consent.

For the identification of metabolomic biomarkers associated with the progression from prediabetes to diabetes, the current study focused on participants with prediabetes at baseline with available circulating metabolite data. The diagnosis of prediabetes was defined by an HbA1c level of 5.7% to 6.4% (39 to 47 mmol/mol) in participants without diabetes, according to the American Diabetes Association (ADA) criteria²². After excluding individuals who developed diabetes or died within 1 month from the baseline, 13,489 participants with prediabetes were included in the final analyses.

Metabolite quantification

The metabolomics analysis of approximately 118,000 non-fasting ethylenediaminetetraacetic acid (EDTA) plasma samples at baseline was performed using the high-throughput NMR platform in Nightingale Health’s laboratories of Finland. Details of the metabolic profiling platform and experimentation have been described elsewhere^23-25. In brief, the EDTA samples were collected and stored at -80°C. Before preparation, frozen samples were slowly thawed at +4°C overnight and were centrifuged (3,400 g) for 3 minutes. Each sample was analyzed with a spectrometer and the metabolic biomarkers were quantified using Nightingale Health’s proprietary software. The quality control procedures were implemented during the whole process and only samples and biomarkers that underwent the quality control process were stored in the UK Biobank dataset.

The metabolic biomarker profiling by Nightingale Health’ s NMR platform provides consistent results over time and across spectrometers. Furthermore, the sample preparation is minimal in the Nightingale Health’ s metabolic biomarker platform, circumventing all extraction steps. These aspects result in highly repeatable biomarker measurements. Pre-specified quality metrics were agreed between UK Biobank and Nightingale Health to ensure consistent results across the samples, and pilot measurements were conducted. Nightingale Health performed real-time monitoring of the measurement consistency within and between spectrometers throughout the UK Biobank samples. Two control samples provided by Nightingale Health were included in each 96-well plate for tracking the consistency across multiple spectrometers. Furthermore, two blind duplicate samples provided by the UK Biobank were included in each well plate, with the position information unlocked only after results delivery. Coefficient of variation (CV) targets across the metabolic biomarker profile were pre-specified for both Nightingale Health’s internal control samples and UK Biobank’s blind duplicates. The targets were met for each consecutively measured batch of ∼25,000 samples. For the majority of the metabolic biomarkers, the CVs were below 5% (https://biobank.ndph.ox.ac.uk/showcase/refer.cgi?id=3000). Further, the distributions of measured biomarkers from 5 sample batches indicated absence of batch effects (https://biobank.ctsu.ox.ac.uk/ukb/ukb/docs/nmrm_app1).

A total of 249 metabolic biomarkers (168 directly measured and 81 ratios of these), spanning lipids, lipoprotein subclass, fatty acids, amino acids, ketone bodies, and glycolysis metabolites were quantified for each sample. In the present study, we analyzed 168 metabolic biomarkers that were directly measured (Supplementary file 1). The values of all metabolites were transformed using natural logarithmic transformation (ln[x+1]) followed by Z-transformation.

Covariate collection

Information on covariates was collected through a self-completed touchscreen questionnaire or verbal interview, including age, sex, ethnicity (white people or others), Townsend deprivation index, household income (high: ≥ £52,000, medium: £18,000-£51,999, and low: < £18,000), education (college/university degree or others), employment status (current working, retired, or other), smoking status (never, previous, or current smoking), moderate alcohol (alcohol intake > 0 g and ≤ 14 g/day for women and alcohol intake > 0 g and ≤ 28 g/day for men), physical activity, healthy diet score, healthy sleep score, family history of diabetes (yes or no), history of cardiovascular disease (CVD, yes or no), history of hypertension (yes or no), history of dyslipidemia (yes or no), history of chronic lung diseases (CLD, e.g. chronic bronchitis, emphysema, and chronic obstructive pulmonary disease, yes or no), and history of cancer (yes or no). The Townsend deprivation index is a composite measure of area-level socioeconomic deprivation, with a higher score indicating higher levels of socioeconomic deprivation. Physical activity was measured by the metabolic equivalent task (MET) (sum of days performing walking, moderate activity, and vigorous activity)²⁷. A healthy diet score was calculated based on the intake of vegetables (≥ median), fruits (≥ median), fish (≥ median), red meat (< median), and processed red meat (< median)²⁸. One point was given for each favorable diet factor and the total diet score ranges from 0 to 5. A healthy sleep score was evaluated based on insomnia (sometimes or never), sleep duration (7∼8 h), chronotype (morning person), daytime sleepiness (sometimes or never), and snoring (no)²⁹. Each favorable sleep factor was given a score of 1, with the total sleep score ranging from 0 to 5. The history of dyslipidemia incorporated information on both the medical history of dyslipidemia and the use of lipid-lowering medications.

Physical measurements including blood pressure, height, weight, waist circumference (WC), and hip circumference (HC) were measured using calibrated instruments with standard protocols by trained nurses. Blood pressure was measured using the Omron automatic digital monitor and two measurements were obtained at a few minutes’ intervals. We calculated the mean systolic (SBP) and diastolic blood pressure (DBP) from two measurements. Body mass index (BMI) was calculated as weight in kilograms divided by the square of height in meters (kg/m²). The HbA1c level was measured by high-performance liquid chromatography with the VARIANT II Turbo analyzer (Bio-Rad Laboratories). Missing covariates were imputed by the median value for continuous variables and a missing indicator for categorical variables.

Ascertainment of diabetes

Incident diabetes was ascertained from hospital inpatient records, death registers, and primary care records, according to the International Classification of Diseases, 10th revision (ICD-10) codes. Detailed information about the linkage procedure is available from https://content.digital.nhs.uk/services. The follow-up time was calculated from the baseline to the occurrence of diabetes, death, or the censoring date (30 March 2023), whichever came first.

Statistical analyses

Baseline characteristics were presented as numbers (percentages) for categorical variables and means (standard deviations, SDs) for continuous variables, respectively. Continuous variables were assessed for statistical differences using t-test and categorical variables were evaluated using the χ² test. Overall schematic workflow of the study is shown in Figure 1.

• Download figure
• Open in new tab

Figure 1. Overall schematic workflow of the study.

CPH, Cox proportional hazard; NMR, nuclear magnetic resonance; RF, random forest; RSF, Random survival forest; SVM, supporting vector machine; XGBoost, extreme gradient boosting.

Metabolite selection

We first used Cox proportional hazards (CPH) model to assess the associations between individual metabolites and risk of diabetes progression with adjustment for sociodemographic covariates (age, sex, ethnicity, education, Townsend Deprivation Index, employment status and household income), family history of diabetes, health conditions (history of CVD, hypertension, dyslipidemia, CLD and cancer), physical measurements (BMI, WC, HC, SBP and DBP), lifestyle factors (smoking status, moderate alcohol, healthy diet score, healthy sleep score and physical activity), and HbA1c. The potential confounders were selected based on prior knowledge of the risk factors for diabetes. Metabolites that were significantly associated with incident diabetes (P <0.05/168) were retained.

Secondly, we performed priority-Lasso to deal with multicollinearity in high dimensional data and to retain variables with nonzero coefficients. Priority-Lasso is a Lasso-based intuitive analysis strategy, which uses prior knowledge regarding the outcome by defining the blocks of different types of predictor variables³⁰. In this study, we defined the 24 covariates as block 1, while all metabolites significantly associated with diabetes risk in the CPH model were defined as block 2. The penalization parameter λ was determined as values with maximum partial-likelihood in a 10-fold cross-validation.

Thirdly, three machine learning models including supporting vector machine (SVM), random forest (RF), and extreme gradient boosting (XGBoost) were adopted to further evaluate the importance of the Lasso-selected metabolites, as they can model nonlinear and nonadditive relations more flexibly³¹. Models were built by 10-fold cross-validation through the “caret” package. Common signals detected across diverse approaches are more likely to represent the strongest and true patterns in the data. We chose the intersection set of the top 20 most important variables selected by the three machine learning models, after balancing the performance of the final diabetes risk prediction model and the clinical applicability associated with measurement costs of metabolites.

Model development

Participants were randomly subclassified into a training set and a test set at a ratio of 8:2 and two common algorithms for survival data including CPH model and random survival forest (RSF)³² were adopted for model development. RSF, as a machine learning method, is designed to be used specifically for survival outcome prediction and has shown promising results in various settings^33,34. It builds many decision trees using split points based on the log-rank test to identify different survival statuses and produces the predicted probability for an individual derived from the average prediction across all trees³⁵. The RSF model was fitted using the “randomForestSRC” package and the grid search method was used for hyperparameter tuning (number of trees, number of variables to possibly split at each node, and minimum size of terminal node). Specifically, the grid search method was used to tune hyperparameters among the RSF model, through minimizing out-of-sample or out-of-bag error³⁶. Each tree in the RSF is constructed from a random sample of the data, typically a bootstrap sample or 63.2% of the sample size (as in the present study). Consequently, not all observations are used to construct each tree. The observations that are not used in the construction of a tree are referred to as out-of-bag observations. In an RSF model, each tree is built from a different sample of the original data, so each observation is “out-of-bag” for some of the trees. The prediction for an observation can then be obtained using only those trees for which the observation was not used for the construction. A classification for each observation is obtained in this way and the error rate can be estimated from these predictions. The resulting error rate is referred to as the out-of-bag error. Through calculating the out-of-bag error in each iteration, the best hyperparameters were finally determined. The hyperparameters to be tuned and range of grid search in the present study were below: number of trees (50-1000, by 50), number of variables to possibly split at each node (3-6, by 1), and minimum size of terminal node (1-20, by 1)³⁷.

Model evaluation

The model performance was assessed in the test set. The time-dependent area under the receiver operating characteristic curve (AUROC) was used to evaluate the model’s discrimination ability. Continuous net reclassification improvement (NRI), and absolute integrated discrimination improvement (IDI) were used to assess whether adding the selected metabolites could improve risk discrimination and reclassification for the risk of progression from prediabetes to diabetes over the basic model that was built on 10 conventional clinical variables (age, sex, Townsend Deprivation Index, family history of diabetes mellitus, BMI, WC, HC, SBP, DBP, and HbA1c)³⁸. The calibration ability of the model was estimated using calibration curve. Furthermore, we used decision curve analysis (DCA) to assess the clinical usefulness of prediction model-based guidance for prediabetes management, which calculates a clinical “net benefit” for one or more prediction models in comparison to default strategies of treating all or no patients³⁹. To facilitate risk stratification, we classified participants into two risk groups according to the predictive value using “surv_cutpoint” function in the “survminer” R package⁴⁰. We also divided participants into three categories according to the tertiles of probability. In addition, we included 90,688 participants with normal glucose from the UK Biobank and divided them into the training and test sets using an 8:2 ratio to further investigate the additive value of the selected metabolites in diabetes prediction among participants with normoglycemia. All analyses were conducted in R software (version 4.2.2). A two-sided P value < 0.05 was considered statistically significant. To control for the false discovery rate in the association between multiple metabolic biomarkers and incident diabetes, Bonferroni correction for P value (P <0.05/168) was used.

Baseline characteristics

Among the 13,489 participants with baseline prediabetes, the mean age was 59.6 (SD, 7.1) years, and 6,166 (45.7%) were males. During a median follow-up of 13.6 (12.3-14.6) years, 2,525 (18.7%) participants progressed to diabetes. Baseline characteristics of the study population stratified by incident diabetes are summarized in Table 1. Participants who developed diabetes were more likely to be male, non-White, less educated, more deprived, and smokers. They also tended to have a family history of diabetes, comorbidities such as CVD, hypertension, dyslipidemia and CLD, and higher levels of BMI, WC, and HC.

Identification of metabolic biomarkers for progression to diabetes

After adjusting for covariates and correcting for multiple testing, 94 of 168 metabolic biomarkers were significantly associated with the risk of incident diabetes (Figure 2, Supplementary file 2). Concentrations of very low-density lipoprotein (VLDL) particles, particularly larger VLDL particles and composition within larger VLDL, were strongly associated with progression to diabetes. Triglyceride in all lipoprotein subclasses also demonstrated strong positive associations with diabetes risk. In contrast, concentrations of larger high-density lipoprotein (HDL) particles and composition within these particles were inversely associated with incident diabetes. For lipoprotein particle diameter, larger HDL and LDL particle sizes were associated with a lower risk of progression to diabetes, while larger VLDL particle size was associated with a higher risk.

• Download figure
• Open in new tab

Figure 2. Associations of 168 metabolic biomarkers with risk of diabetes among 13,489 participants with prediabetes.

Hazard ratios (HR) were presented per 1 standard deviation (SD) higher of metabolic biomarker on the natural log scale and were adjusted for age, sex, ethnicity, education, Townsend Deprivation Index, employment status, household income, family history of diabetes, history of CVD, history of hypertension, history of dyslipidemia, history of CLD, history of cancer, body mass index, waist circumference, hip circumference, smoking status, moderate alcohol, healthy diet score, healthy sleep score, physical activity, systolic blood pressure, diastolic blood pressure and glycated hemoglobin A1c. *False discovery rate controlled P < 0.05/168.

Apo-A1, apolipoprotein A1; Apo-B, apolipoprotein B; Apo-LP, apolipoprotein; BCAA, branched-chain amino acid; BMI, body mass index; CVD, cardiovascular disease; CLD, chronic lung disease; DHA, docosahexaenoic acid; FA, fatty acids; HDL, high-density lipoproteins; HDL-D, high-density lipoprotein particle diameter; IDL, intermediate-density lipoproteins; L, large; LA, linoleic acid; LDL, low-density lipoproteins; LDL-D, low-density lipoprotein particle diameter; LP, lipoprotein; M, medium; MUFA, monounsaturated fatty acids; PUFA, polyunsaturated fatty acids; S, small; SFA, saturated fatty acids; VLDL, very-low-density lipoproteins; VLDL-D, very-low-density lipoprotein particle diameter; XL, very large; XS, very small; XXL, extremely large.

Monounsaturated fatty acids and saturated fatty acids were positively associated with the risk of diabetes, whereas docosahexaenoic acid and the degree of fatty acid unsaturation were negatively associated with diabetes. Among the amino acids, higher concentrations of alanine, tyrosine, and branched-chain amino acid (BCAA) such as leucine and valine were associated with an increased risk of diabetes, but glutamine and glycine were inversely associated with diabetes. Neither of the ketone bodies showed an association with the risk of diabetes.

Of the 94 metabolites that were significantly associated with diabetes, 17 metabolites were selected by priority-Lasso (Supplementary file 3). When further evaluating the importance of these metabolites after adjustment for covariates using three machine learning algorithms, the intersection of the top 20 important predictors identified a total of 9 metabolites, namely cholesteryl esters in large HDL, cholesteryl esters in medium VLDL, triglycerides in very large VLDL, average diameter for LDL particles, triglycerides in IDL, glycine, tyrosine, glucose, and docosahexaenoic acid (Figure 3, Supplementary file 4).

• Download figure
• Open in new tab

Figure 3. The top 20 important variables selected by three machine learning models: (A) supporting vector machine (SVM); (B) extreme gradient boosting (XGBoost); (C) random forest (RF).

The models were adjusted for age, sex, ethnicity, education, Townsend Deprivation Index, employment status, household income, family history of diabetes, history of CVD, history of hypertension, history of dyslipidemia, history of CLD, history of cancer, body mass index, waist circumference, hip circumference, smoking status, moderate alcohol, healthy diet score, healthy sleep score, physical activity, systolic blood pressure, diastolic blood pressure and glycated hemoglobin A1c. CVD, cardiovascular disease; CLD, chronic lung disease. HDL, high-density lipoproteins; IDL, intermediate-density lipoproteins; LDL, low-density lipoproteins; VLDL, very-low-density lipoproteins.

Model development and evaluation

Build upon the selected 9 metabolites and 10 clinical variables, there was no obvious difference in the AUROC obtained from CPH model (1-year: 0.823 [95% confidence interval, CI 0.702, 0.945]; 5-year: 0.830 [0.797, 0.864]; 10-year: 0.801 [0.778, 0.825]) and RSF model (1-year: 0.828 [0.723, 0.933]; 5-year: 0.820 [0.785, 0.855]; 10-year: 0.802 [0.778, 0.826]). Hence, we chose CPH model as the final model because of its simplicity and interpretability. The addition of selected metabolites consecutively outperformed the basic model with conventional clinical variables in diabetes risk prediction from 1 to 10 years (Figure 4). Specifically, the AUROC increased from 0.759 (95% CI 0.608, 0.911) to 0.823 (0.702, 0.945), 0.798 (0.762, 0.834) to 0.830 (0.797, 0.864), and 0.776 (0.750, 0.801) to 0.801 (0.778, 0.825) for 1-year, 5-year, and 10-year diabetes risk, respectively (Table 2, Figure 5). Results from continuous NRI and absolute IDI also demonstrated improvement in the risk prediction for progression to diabetes (Table 2), although the model calibration was not significantly improved (Figure 6). The decision curve analysis showed that the inclusion of the metabolites had a higher net benefit across the threshold probabilities of 0-0.35 for predicting 5-year diabetes risk and 0-0.55 for predicting 10-year diabetes risk (Figure 7).

• Download figure
• Open in new tab

Figure 4. Consecutive area under time-dependent receiver-operating characteristic (AUROC) of basic model and basic model plus 9 metabolites (A), and the difference of these two time-dependent AUROCs over time (B).

The basic model used conventional clinical variables including age, sex, Townsend Deprivation Index, family history of diabetes mellitus, body mass index, Waist circumference, hip circumference, systolic blood pressure, diastolic blood pressure, and glycated hemoglobin A1c. The selected 9 metabolites included cholesteryl esters in large HDL, triglycerides in very large VLDL, Glycine, average diameter for LDL particles, tyrosine, cholesteryl esters in medium VLDL, glucose, triglycerides in IDL, and docosahexaenoic acid.

HDL, high-density lipoprotein; IDL, intermediate-density lipoprotein; LDL, low-density lipoprotein; VLDL, very-low-density lipoprotein.

• Download figure
• Open in new tab

Figure 5. Time-dependent receiver-operating characteristic (ROC) curves of basic model and basic model plus 9 metabolites for predicting 1-year (A), 5-year (B), and 10-year (C) risk of developing diabetes in participants with prediabetes.

The basic model used conventional clinical variables including age, sex, Townsend Deprivation Index, family history of diabetes mellitus, body mass index, Waist circumference, hip circumference, systolic blood pressure, diastolic blood pressure, and glycated hemoglobin A1c. The selected 9 metabolites included cholesteryl esters in large HDL, triglycerides in very large VLDL, Glycine, average diameter for LDL particles, tyrosine, cholesteryl esters in medium VLDL, glucose, triglycerides in IDL, and docosahexaenoic acid.

HDL, high-density lipoprotein; IDL, intermediate-density lipoprotein; LDL, low-density lipoprotein; VLDL, very-low-density lipoprotein.

• Download figure
• Open in new tab

Figure 6. Calibration plots of basic model (A, B, C) and basic model plus 9 metabolites (D, E, F) for predicting 1-year, 5-year, and 10-year risk of developing diabetes in participants with prediabetes.

The basic model used conventional clinical variables including age, sex, Townsend Deprivation Index, family history of diabetes mellitus, body mass index, Waist circumference, hip circumference, systolic blood pressure, diastolic blood pressure, and glycated hemoglobin A1c. The selected 9 metabolites included cholesteryl esters in large HDL, triglycerides in very large VLDL, Glycine, average diameter for LDL particles, tyrosine, cholesteryl esters in medium VLDL, glucose, triglycerides in IDL, and docosahexaenoic acid.

HDL, high-density lipoprotein; IDL, intermediate-density lipoprotein; LDL, low-density lipoprotein; VLDL, very-low-density lipoprotein.

• Download figure
• Open in new tab

Figure 7. Decision curve analysis of basic model and basic model plus 9 metabolites for predicting 5-year (A) and 10-year (B) risk of developing diabetes in participants with prediabetes.

Decision curve analysis was not performed on 1-year prediction considering the relatively small number of prediabetic patients who develop diabetes within a year in the test set and small net benefit from intervention. The basic model used conventional clinical variables including age, sex, Townsend Deprivation Index, family history of diabetes mellitus, body mass index, Waist circumference, hip circumference, systolic blood pressure, diastolic blood pressure, and glycated hemoglobin A1c. The selected 9 metabolites included cholesteryl esters in large HDL, triglycerides in very large VLDL, Glycine, average diameter for LDL particles, tyrosine, cholesteryl esters in medium VLDL, glucose, triglycerides in IDL, and docosahexaenoic acid.

HDL, high-density lipoprotein; IDL, intermediate-density lipoprotein; LDL, low-density lipoprotein; VLDL, very-low-density lipoprotein.

We further categorized the participants from the test set into low-risk and high-risk groups according to the optimal threshold of the predicted value (1.02) reflecting the best risk difference. Compared with the low-risk group, participants in the high-risk group had a significantly higher cumulative risk of incident diabetes (log-rank P < 0.0001) (Figure 8). When participants were alternatively classified into low-risk, medium-risk, and high-risk groups according to the tertile cut-off point of the predicted value, the high-risk group showed the highest risk of developing diabetes, followed by the medium-risk and low-risk groups (log-rank P < 0.0001). Similar results were also observed when considering the competing risk from death (Fine-Gray P < 0.0001) (Figure 8-figure supplement 1). In addition, the predicted risk of diabetes within 1 year (P = 0.001), 5 years (P < 0.001), or 10 years (P < 0.001) was generally higher among participants who progressed to diabetes than those who did not (Figure 9).

• Download figure
• Open in new tab

Figure 8-figure supplement 1. Cumulative hazard curves for participants with prediabetes with different risks stratified by the Cox model based on clinical variables and 9 metabolites when considering competing risk from death.

The Cox model divided participants with prediabetes in the test set to two categories (A) and three categories (B) with significant differences in cumulative hazard of diabetes during the follow-up (both Fine-Gray P <0.0001).

• Download figure
• Open in new tab

Figure 8. Cumulative hazard curves for participants with prediabetes with different risks stratified by the Cox model based on clinical variables and 9 metabolites.

The Cox model divided participants with prediabetes in the test set to two categories (A) and three categories (B) with significant differences in cumulative hazard of diabetes during the follow-up (both P <0.0001).

• Download figure
• Open in new tab

Figure 9. The distribution of the predictive probability of developing diabetes among participants with prediabetes by incident diabetes status within 1-year (A), 5-year (B), and 10-year (C).

Among participants with normoglycemia, we also observed a significant improvement in the prediction of diabetes after the addition of metabolic biomarkers to the basic model. The AUROC increased from 0.821 (95% CI 0.736, 0.907) to 0.868 (0.802, 0.934), 0.790 (0.738, 0.842) to 0.811 (0.762, 0.860), and 0.791 (0.765, 0.816) to 0.806 (0.781, 0.831) for 1-year, 5-year, and 10-year diabetes risk, respectively (Supplementary file 5). The increases in NRI and IDI were similar to or slightly lower than those found among participants with prediabetes.