J. Diana Zhang,Chonghua Xue,Vijaya B. Kolachalama,William A. Donald
Interpretable Machine Learning on Metabolomics Data Reveals Biomarkers for Parkinson’s Disease
The use of machine learning (ML) with metabolomics provides opportunities for the early diagnosis of disease. However, the accuracy of ML and extent of information obtained from metabolomics can be limited owing to challenges associated with interpreting disease prediction models and analyzing many chemical features with abundances that are correlated and “noisy”. Here, we report an interpretable neural network (NN) framework to accurately predict disease and identify significant biomarkers using whole metabolomics data sets without a priori feature selection. The performance of the NN approach for predicting Parkinson’s disease (PD) from blood plasma metabolomics data is significantly higher than other ML methods with a mean area under the curve of >0.995. PD-specific markers that predate clinical PD diagnosis and contribute significantly to early disease prediction were identified including an exogenous polyfluoroalkyl substance. It is anticipated that this accurate and interpretable NN-based approach can improve diagnostic performance for many diseases using metabolomics and other untargeted ‘omics methods.
Interpretable Machine Learning on Metabolomics Data Reveals Biomarkers for Parkinson’s Disease

The use of machine learning (ML) with metabolomics provides opportunities for the early diagnosis of disease. However, the accuracy of ML and extent of information obtained from metabolomics can be limited owing to challenges associated with interpreting disease prediction models and analyzing many chemical features with abundances that are correlated and “noisy”. Here, we report an interpretable neural network (NN) framework to accurately predict disease and identify significant biomarkers using whole metabolomics data sets without a priori feature selection. The performance of the NN approach for predicting Parkinson’s disease (PD) from blood plasma metabolomics data is significantly higher than other ML methods with a mean area under the curve of >0.995. PD-specific markers that predate clinical PD diagnosis and contribute significantly to early disease prediction were identified including an exogenous polyfluoroalkyl substance. It is anticipated that this accurate and interpretable NN-based approach can improve diagnostic performance for many diseases using metabolomics and other untargeted ‘omics methods.

The rate of Parkinson’s disease (PD) is growing more rapidly than any other neurological disease. PD is typically diagnosed according to a clinical criteria of motor symptoms which include bradykinesia (slowness of movement), a resting tremor, and rigidity. However, the onset of atypical nonmotor symptoms such as sleep disorder, constipation, apathy, and loss of smell can predate clinically relevant symptoms by several years to decades. In addition, for patients who present with Parkinson-like symptoms, the current process for identifying PD can often be inconclusive. For example, according to a meta-analysis by Rizzo et al., the overall diagnostic accuracy for PD based on an initial clinical assessment by movement disorder experts is 80%. Accurate identification of PD using biomarker signatures rather than relying primarily on physical symptoms would be highly beneficial. Biomarkers associated with metabolic processes are used extensively for understanding, diagnosing, and monitoring diseases. Such metabolites are typically sampled from well-established matrices such as blood plasma and serum for trace-level analysis of up to thousands of metabolites using mass spectrometry (MS). Additional matrices of emerging interest for biomarker discovery and disease diagnosis applications include the rapid and noninvasive sampling of skin sebum and breath. Using MS, differences in the metabolite profiles in the blood plasma of pre-PD subjects were identified up to 15 years prior to a clinical diagnosis when compared to healthy controls who did not develop PD. These results suggest that PD may potentially be diagnosed using metabolite biomarkers significantly earlier than in current practice, particularly if analyzing such metabolites can result in high diagnostic accuracy and be validated on large-scale cohort studies. To develop accurate prediction models for disease diagnosis using large metabolomics data sets, machine learning (ML) approaches are widely used. However, the use of whole metabolomics data sets to build prediction models is rare. Such data sets can contain highly correlated and “noisy” chemical features which may risk model overtraining and reduce diagnostic performance. As a result, models are typically based on a smaller subset of features which are determined by conventional statistical methods (e.g., based on p-values and fold-changes of individual features). For example, Gonzalez-Riano et al. used a linear support vector machine (SVM) model with 20 preselected biomarkers to diagnose pre-PD vs healthy controls from blood plasma samples. Similarly, Sinclair et al. used partial least-squares-discriminant analysis (PLS-DA) with 15 and 26 preselected biomarkers to diagnose drug naïve PD and medicated PD vs healthy control, respectively, from skin sebum samples. However, given that the abundances of metabolites are often correlated and can depend nonlinearly on the abundances of other metabolites, ML approaches such as SVM and PLS-DA may potentially “miss” some key features in metabolomics data sets. Advanced ML approaches such as neural networks (NN) are particularly well-suited for processing large volumes of correlated data and building models for data sets that contain nonlinear effects. However, a fundamental issue in using methods such as NN for classifying complex mixtures based on metabolomics data is that the resulting predictive models are generally considered as uninterpretable “black boxes”, which cannot be readily used to reveal mechanistic information. Recently, a new approach entitled Shapley Additive exPlanations (SHAP) was developed to “interpret” ML models by retrospectively calculating the contribution of individual features to the accurate predictive performance of a model. However, SHAP has not been used in the analysis of metabolomics data sets given that methods for interpreting ML models have only recently been developed, and using all chemical features risks overtraining prediction models. Ideally, whole metabolomics data sets should be included in the ML model for SHAP to identify key metabolites that drive model prediction. Here, we report an interpretable neural network-based framework for analyzing data sets generated by untargeted mass spectrometry-based methods (Figure 1) entitled, “CRANK-MS” (Classification and Ranking Analysis using Neural network generates Knowledge from Mass Spectrometry). CRANK-MS has several built-in features including (i) integrated model parameters that allow the high dimensionality of metabolomics data sets to be analyzed without the need for preselecting chemical features; (ii) SHAP to retrospectively “mine” key chemical features that contribute the most to an accurate model prediction; and (iii) benchmark testing with five well-known ML methods to compare diagnostic performance and further verify significant chemical features. Using CRANK-MS, we report the highest diagnostic performance to date for binary classification of PD vs healthy control. Additionally, NN-driven predictions trained on a prognostic PD study were used to reveal new PD-specific chemical features which were not previously identified and can be considered indicative of pre-PD diagnosis. The program for implementing this approach is freely available online at https://github.com/CRANK-MS.

Data sets from two cross-sectional PD metabolomics studies were used throughout. The Spanish European Prospective Study on Nutrition and Cancer (EPIC) study involved metabolomics data from blood plasma samples taken from subjects who later developed PD up to 15 years later, and those who did not develop PD (total number of participants, n = 78; Table 1). The blood plasma samples from the EPIC study were analyzed using four different instrumental methods (gas chromatography-MS, GC-MS; capillary electrophoresis-MS, CE-MS; and liquid chromatography-MS, LC-MS, in positive (+) and negative (−) ionization modes). The NHS study involved the LC-MS (+) analysis of skin sebum sampled from drug-naïve and medicated PD patients, and healthy controls (n = 274). The composite data set from the EPIC study was prepared using the data from all four methods (i.e., GC-MS, CE-MS, LC-MS (+), and LC-MS (−)). Six participants were excluded in the composite data set as data from one or more of the four methods were missing. The number of reported molecular features in the metabolomics data sets ranged from 60 to 6502 (Table 1).

Using the metabolomics data as inputs, we leveraged six supervised learning frameworks to classify persons with PD from those who are healthy. Random forest (RF), extreme gradient boosting (XGB), linear discriminant analysis (LDA), logistic regression (LR), and support vector machine (SVM) were written using scikit-learn packages (v 1.0.2). The algorithms and SHAP analysis were implemented in Python (v. 3.8). NN with multilayer perceptron was written using PyTorch (v. 1.10.2). Additional Python libraries used to support data analysis and visualization include pandas (v. 1.4.2), numpy (v. 1.21.5), and matplotlib (v. 3.5.1). Full details of the open access code are available at https://github.com/CRANK-MS.

The pseudocode for each ML algorithm based on 100 bootstraps is presented in Chart 1.

For each ML model, hyperparameters correspond to specific model parameters that an algorithm uses to train on a given data set. To determine the optimal hyperparameters, the composite data set from the EPIC study was used (see above). Hyperparameter tuning for each ML model was optimized using the GridSearchCV package in scikit-learn. For NN, XGB, RF, LR, and SVM, the number of possible permutations resulting from different parameters ranged from 120 to 158. For LDA, the number of permutations used was 12. For each permutation, a bootstrap model was used in which the data set was split randomly 100 times into 60% training data and 40% validation data (i.e., 100 “bootstraps”). The optimized hyperparameters were determined based on the combination of hyperparameters that resulted in the highest Matthews correlation coefficient (MCC) after 100 iterations per permutation. These hyperparameters were then applied to all data sets in the study (Table S1).

To calculate diagnostic performance, the 100 randomly selected testing data sets from the bootstrapping process were used to calculate the overall diagnostic performance. The use of a bootstrap model with more replicates can lower absolute error compared to other sampling methods such as cross-validation and is considered useful for relatively small sample sizes. The diagnostic performance for each ML model was calculated based on the mean of the 100 bootstrap measurements, and error was calculated as one standard deviation of the mean. For each ML model, accuracy, precision, sensitivity/recall, specificity, F1 score, and MCC score were calculated. Receiver operating characteristic (ROC) and precision-recall (PR) curves were generated to calculate area-under-curve (AUC). Briefly, AUC (ROC) is a plot of the true positive rate (i.e., how many PD patients were correctly predicted) vs the true negative rate (i.e., how many healthy controls were correctly predicted). In contrast, AUC (PR) plots the precision rate (i.e., how many PD predictions were correct) vs recall or sensitivity rate (i.e., how many PD patients were correctly predicted).

For each chemical feature in the metabolomics data sets, a SHAP score was calculated based on the absolute average from the 100 bootstraps. The greater the SHAP scores, the greater the contribution to the overall prediction of PD across all 100 bootstraps and all participants. Metabolites were annotated based on the highest mass accuracies that were obtained by comparing the measured monoisotopic neutral masses (accounting for protonation, sodiation, and potential loss of a water molecule) to those from the Human Metabolome Database (https://hmdb.ca/) using a threshold of ±20 ppm. One top scoring chemical feature (m/z 942.9824) had a relatively large negative mass defect which is indicative of an exogenous, synthetic compound. Thus, this ion was annotated using the PubChem (http://www.cheminfo.org/) database.

The overall diagnostic performance for all six ML algorithms was assessed using a composite data set featuring metabolites from blood plasma that were detected using four analytical methods as reported in the EPIC study. The diagnostic performance of NN was higher than the other frameworks across all metrics investigated (Figure 2, Table S2). Specifically, the binary classification of PD vs healthy using NN resulted in AUCs of 0.994 ± 0.018 and 0.995 ± 0.014 for ROC and PR, respectively. Extreme gradient boosting and logistic regression performed similarly with AUC (ROC) and AUC (PR) of 0.970 ± 0.028 and 0.968 ± 0.031 for extreme gradient boosting, and 0.968 ± 0.037 and 0.969 ± 0.037 for logistic regression, respectively. In contrast, the performance of the RF, SVM, and LDA classifiers were relatively low with AUC (ROC) and AUC (PR) values of 0.829 ± 0.099 and 0.836 ± 0.099 for the RF classifier, 0.647 ± 0.093 and 0.661 ± 0.111 for the SVM classifier, and 0.681 ± 0.091 and 0.634 ± 0.119 for the LDA classifier, respectively. Additional metrics, i.e., accuracy, sensitivity, specificity, precision, and F1 score can be found in Table S2. For all six algorithms, randomly permuting the PD and healthy data labels resulted in average accuracies (Figure S1) that are statistically the same as 0.5 (50%), consistent with a “random guess” for binary classification as expected for this control test. Recently, Chicco et al. has shown that the MCC score is a more informative and reliable metric for evaluating binary classification accuracy as it considers all four values in the confusion matrix (i.e., true positive, false positive, true negative, and false negative). Thus, MCC can be considered as a less biased metric toward data sets with imbalanced cohorts. Based on the MCC score, NN performed significantly higher with 0.918 ± 0.086 compared to 0.815 ± 0.132, 0.787 ± 0.119, 0.433 ± 0.192, 0.272 ± 0.152, and 0.213 ± 0.155 for the LR, XGB, RF, LDA, and SVM classifiers, respectively. SHAP analysis was used to identify the metabolites and the corresponding mass spectrometry-based method that contributed the most to the prediction of PD using the composite metabolomics data set for blood plasma. Five of the top six metabolites were detected using LC-MS (+) (i.e., m/z 942.9824, 467.3822, 393.3454, 379.3289, and 613.4767) (Figure 2). To further validate the contribution of these chemical features in predicting PD, all six ML algorithms were applied to the LC-MS (+) data set without including the LC-MS (−), GC-MS, and CE-MS data sets. Similar to that for the composite data set, binary classification of PD using NN and the LC-MS (+) data was highest across all performance metrics (Figures 3 and S2, Table S2). For example, the AUC (ROC), AUC (PR), and MCC score for NN was 0.983 ± 0.022, 0.984 ± 0.022, and 0.894 ± 0.081, respectively. Both XGB and LR classifiers performed about the same or slightly lower with AUC (ROC), AUC (PR), and MCC values of 0.968 ± 0.029, 0.966 ± 0.036, 0.805 ± 0.117, and 0.972 ± 0.029, 0.976 ± 0.026, 0.869 ± 0.078, respectively. RF, SVM, and LDA classifiers were substantially lower with AUC (ROC), AUC (PR), and MCC scores of 0.894 ± 0.079, 0.911 ± 0.071, 0.589 ± 0.177 for the RF classifier; 0.856 ± 0.073, 0.869 ± 0.078, 0.582 ± 0.167 for the SVM classifier; and 0.878 ± 0.066, 0.877 ± 0.086, and 0.626 ± 0.150 for the LDA classifier. Additional metrics are given in Table S2. The standard deviations for all performance metrics were pooled to calculate the average relative standard deviation (RSD). The RSD was lowest for NN (5.17%) followed by LR (5.19%), XGB (7.70%), LDA (12.8%), SVM (14.3%), and RF (17.2%). Based on a SHAP analysis, five of the six top-scoring chemical features in the LC-MS (+) data set were also in the top six highest scoring features for the composite data set (see above, Figure 3). These results indicate that the diagnostic performance for the LC-MS (+) data set outperformed the three other mass spectrometry-based methods across all six ML methods. The use of NN resulted in higher diagnostic accuracy than the other five ML methods for the LC-MS (+) data set, which has significantly fewer chemical features (509) than in the composite data set (1430). Overall, NN resulted in the highest diagnostic performance and the lowest RSD in predicting PD from blood plasma using either the composite or LC-MS (+) data sets compared to the other five ML algorithms. Disease classification using NN involving 100 bootstraps and 1430 total metabolites required <1 min on a consumer laptop computer (Surface Laptop 3, Microsoft) with a 1.2 GHz processor (Intel i5-core). The diagnostic performance of NN was at least 10% higher than previously reported by Gonzalez-Riano et al. when a SVM model was used. In addition, the diagnostic performance obtained using NN is the highest reported to date for any PD diagnosis regardless of the sample matrix including blood plasma, blood serum, and skin sebum.

The NN model included all metabolites or features from the data set as inputs unlike training models with preselected chemical features. For example, in Gonzalez-Riano et al., biomarkers for PD were first screened for significance, and a small subset of these biomarkers (up to 20) was used in the final diagnostic model. In this previous study, the highest ROC (AUC) value obtained was 0.919 for the composite data set using a 20-feature linear SVM model. In the current study, a similar model was applied to the composite data set without feature selection which resulted in an AUC (ROC) of 0.647 ± 0.093. In contrast, using NN on all features in the composite data set resulted in an AUC (ROC) of 0.994 ± 0.018. These results are consistent with some well-known ML models having relatively low predictive performance when incorporating large data sets that contain many “noisy” features. A feature-selected NN model was developed using the ten chemical features that had the highest SHAP scores in the LC-MS (+) data set. The AUC ROC and PR values for the feature-selected NN-model were comparable or slightly higher (0.997 ± 0.006 and 0.997 ± 0.006) than that obtained using all chemical features in the data set (0.983 ± 0.022 and 0.984 ± 0.022). Given that the diagnostic performance of both models was comparable, these data indicate that NN can be highly tolerant of many chemical features (>1500) that do not contribute substantially to accurate disease prediction; i.e., diagnostic performance is essentially uncompromised by including the whole metabolomics data set without feature selection. In addition, the relatively high diagnostic accuracy supports the use of SHAP to accurately identify chemical features that contribute significantly to disease classification.

The metabolites that contributed the most significantly to the accurate prediction of PD are more likely to be more basic and readily ionized by cation adduction, rather than acidic, given that higher diagnostic performance for PD was obtained using LC-MS (+) than LC-MS (−) and that similar numbers of metabolites (∼510 to 530) were measured using each method, consistent with previous reports. SHAP analysis on LC-MS (+) data for blood plasma revealed that five of the top six highest scoring metabolites were consistent across all six ML algorithms (Figure 3). The detected metabolites were different compared to those determined using a previous study focused on a linear SVM model (Figure S3), which can be attributed to the difference between using kernel-based methods and neural networks. The five metabolites that contributed the most to an accurate prognostic PD prediction could serve as potential indicators for disease status and were annotated (Table 2). The five annotated ions corresponded to a polyfluoroalkyl substance (PFAS), triterpenoids, cholestane steroids, diacylglycerol, and vitamin D steroids of either endogenous or exogenous origins, which have been linked to PD in the literature previously (Table 2). For example, the ion with an m/z value of 942.9824 had the highest SHAP value and likely corresponds to the sodiated PFAS [3-(2,2,3,3,4,4,4-heptafluorobutanoyloxy)-2,2-bis(2,2,3,3,4,4,4-heptafluorobutanoyloxymethyl)propyl] 2,2,3,3,4,4,4-heptafluorobutanoate (DTXSID70325550). Ions corresponding to this PFAS and its oxidation ([3-(2,2,3,3,4,4,4-heptafluorobutanoyloxy)-2,2-bis(2,2,3,3,4,4,4-heptafluorobutanoyloxymethyl)-1-hydroxypropyl] 2,2,3,3,4,4,4-heptafluorobutanoate) and hydrolysis (i.e., 2,2-bis(hydroxymethyl)propane-1,3-diol) products were all higher in PD participants than healthy controls. The presence of PFAS compounds are ubiquitous in the environment and human blood given their propensity to bioaccumulate, chemical longevity, and widespread use in industrial and consumer products such as plastics, nonstick cookware, and food packaging. For example, in the U.S. population, PFAS was detected in the blood serum of over 98% of Americans that were sampled during 2003–2004 (n = 2,094). DTXSID70325550 is a PFAS compound of interest that is currently listed under the U.S. Environmental Protection Agency CompTox Chemicals Database and appears preorganized for the noncovalent, multidentate binding of Na+, K+, Ca2+, Cu2+, and Zn2+. Thus, such a compound could potentially disrupt neuronal activity by affecting intracellular ion homeostasis. A potential mechanism proposed for PFAS-induced neurotoxicity involves the increase of intracellular Ca2+ which is implicated in impacting neuronal cell processing, signaling, and function. Although further in vitro and in vivo studies are needed to investigate the effects of DTXSID70325550 on neuronal cell function, these data suggest that elevated levels of specific PFAS compounds in blood plasma may be an early indicator of PD. Overall, these results further support that SHAP analysis can be useful in identifying potential biomarkers for PD (Table 2) that were not initially found using statistical approaches.

The performance of CRANK-MS was assessed on a larger sample size with more chemical features. The six ML algorithms were used to analyze the data from the NHS study to predict PD patients that were drug-naïve or medicated from healthy controls using LC-MS (+) metabolomics data from skin swab samples of sebum. On average across all performance metrics, the performance of NN was higher than the five other algorithms (Figure 4, Table S2). Binary classification of drug-naïve PD vs healthy control using NN resulted in AUC (ROC), AUC (PR), and MCC scores of 0.843 ± 0.045, 0.896 ± 0.046, and 0.530 ± 0.098 respectively. Additional performance metrics are shown in Table S2. In addition, the average RSD was lowest for NN at 9.51% followed by SVM (9.63%), LR and LDA (10.0%), XGB (11.7%), and RF classifiers (14.0%). The accuracy in predicting medicated PD vs healthy control across all six algorithms was about the same or slightly higher compared to drug-naïve PD vs healthy control (Table S2), consistent with the medicated PD cohort being more accurately diagnosed based on a clinical assessment of their symptoms and/or medication influencing skin sebum metabolites. For the drug-naïve PD vs healthy control data set, the diagnostic performance obtained using NN on skin sebum was more than 8% higher than previously reported by Sinclair et al. which used a multivariate principle least-squares-discriminant analysis model based on 15 preselected features. Specifically, the AUC (ROC) using multivariate principle least-squares-discriminant analysis was 0.779 compared to 0.843 ± 0.045 using NN. The NN approach required no preselection of features and took 10 min to obtain the key diagnostic performance metrics (i.e., accuracy, sensitivity, specificity, precision, F1 score, MCC score, AUC ROC, AUC PR) using all 6,502 chemical features in the data set. The relatively lower diagnostic performance obtained for the NHS study compared to the EPIC study may be attributable to different study designs (Table 1) in addition to differences in the use of blood plasma vs skin sebum. For example, in the EPIC study, PD participants were matched to healthy subjects based on several factors including the same recruitment center, age, and time of day for sample collection, and fasting period, unlike in the NHS study.

In this study, a NN framework with interpretable feature analysis, entitled CRANK-MS, is reported that can be used to establish accurate disease prediction models using whole MS data sets without preselecting features. The NN is highly tolerant of “noisy” metabolomics data that can contain thousands of metabolites which may not contribute significantly to model prediction. Using CRANK-MS, we report the highest diagnostic performance to date for predicting PD using blood plasma metabolomics data (0.997 AUC for both ROC and PR) when benchmarked with well-known ML algorithms such as PLS-DA and SVM. Diagnostic accuracy in predicting PD using skin sebum metabolomics data was also enhanced using NN compared to alternative, widely used ML approaches. PD-specific biomarkers including triterpenoids, diacylglycerols, and a polyfluoroalkyl substance that contributed significantly to ML model predictions were identified from blood samples that were collected up to 15 years prior to when subjects were clinically diagnosed with PD. These data indicate that these metabolites are potential early indicators for PD that predate clinical PD diagnosis and are consistent with specific food diets (such as the Mediterranean diet) for PD prevention and that exposure to some exogenous chemicals (such as PFASs that can disrupt neuronal activity via changes to intracellular ion homeostasis) may contribute to the development of PD. Given the improved diagnostic performance of CRANK-MS, it is anticipated that this NN-based framework can be a powerful tool to build accurate prediction models for other diseases using metabolomics data. Interpretable ML methods can also be used to retrospectively “mine” metabolomics data sets to identify early “lead” compounds within the biomarker discovery pipeline. Biomarkers identified using this approach can be further validated using high-resolution tandem mass spectrometry and NMR for complete structure elucidation and targeted quantification using clinical MS-based methods, in addition to using in vitro cell-based assays and in vivo disease models. For example, future studies on targeted quantification of pre-PD metabolites could be used to track and monitor disease progression by integrating data on clinical deterioration and other factors such as lifestyle, medication, and diet. In addition, given that there are over 800 publicly available metabolomics studies in The Metabolomics Workbench data repository (https://www.metabolomicsworkbench.org/), prediction models could be used for the binary classification of diseases such as diabetes, fatty liver disease, heart disease, chronic obstructive pulmonary disease, and COVID-19. Using advanced ML methods, retrospectively “mining” such databases for biomarkers that contribute significantly to the prediction of these diseases could reveal novel mechanistic information that may not necessarily be apparent using traditional linear approaches. In addition, CRANK-MS could be adapted to support clinical workflows and improve confidence in diagnosis including for diseases in which stratification is important by performing multinomial classification (>2 classes). For example, CRANK-MS can be used with metabolomics data in conjunction with alternative clinical information such as medical history and neurological examination scores, as well as neuroimaging, to further differentiate clinical PD from other types of Parkinsonian-like diseases. Ultimately, the use of CRANK-MS should enhance the accuracy of disease prediction models based on metabolomics and many other types of ‘omics experiments and facilitate biomarker discovery.