Harnessing Metabolites as Serum Biomarkers for Liver Graft Pathology Prediction Using Machine Learning

1. Introduction

Chronic graft injury compromises long-term survival: Liver transplantation saves thousands of lives worldwide annually – in fact, there were over 37,000 transplants in 2022, and 13,400 in the Americas alone[1]. However, long-term survival in 25% of liver transplant (LT) recipients is compromised by ongoing graft injury that results in cirrhosis[2]. Liver graft injury is typically heralded by abnormalities in liver biochemistry[3,4]. Causes of graft injury include T-cell mediated rejection (TCMR), metabolic dysfunction-associated steatohepatitis (MASH), biliary complications and viral infections, amongst others[5]. TCMR is the most common cause of liver graft injury, and repeated episodes lead to chronic rejection, premature graft loss and compromised long-term survival[6]. MASH recurs in most patients originally transplanted for MASH (MASH-LT), and an estimated 50% develop significant graft fibrosis (defined as Stage 2 or greater) within 5 years of transplant[7]. Overall, ongoing graft injury leads to accelerated fibrosis in comparison to the native liver, progressing at an estimated rate of 0.4 stages per year, which can rapidly result in cirrhosis and loss of the graft[6,8,9,10,11].

How can we best preserve the long-term health of the liver graft? It is imperative that we optimize the long-term outcomes of the graft and its recipient using a personalized, data-driven approach. The only way to reliably diagnose graft injury at present is by performing a liver biopsy and assessing histological features. However, a liver biopsy is an invasive procedure with a 1.8% risk of complications, and it is impractical to perform longitudinal liver biopsies over a LT recipient’s lifetime[12]. Thus, there is great need for effective serological biomarkers to facilitate the noninvasive diagnosis of rejection, biliary complications and MASH.

The term ‘metabolomics’ describes the identification and quantification of metabolites in biological tissue[13]. Metabolites, as downstream products of gene expression, protein and enzymatic function can provide valuable information on the biological processes within a cell, tissue, organ system or organism in addition to the pathophysiology behind different disease states[13,14]. In the pre-transplant population, metabolomics studies have identified potential non-invasive biomarkers associated with advanced MASH and hepatocellular carcinoma[15,16,17,18]. However, data is quite limited, especially in post-transplant populations, although early single center studies have identified potential metabolomic profiles in donors associated with early allograft dysfunction, indicators of ischemia-reperfusion injury in recipients and TCMR in pediatric populations[14,19,20].

This study aimed to identify distinct metabolomic profiles in the serum of individuals following liver transplant, with the goal of recognizing potential biomarkers capable of differentiating between post-transplant complications, specifically MASH, TCMR, and biliary issues. Serum metabolomics is relatively cost-effective as compared to other high-throughput approaches. We employed a random forest (RF) algorithm to develop a classification model distinguishing between biliary complications, MASH, and TCMR. The interpretation of the model was conducted through the permutation-based feature importance measurement for random forests[21].

2. Materials and Methods

2.1. Sample Collection and Processing

Serum samples from consented patients with MASH (n=10), TCMR (n=18), biliary complications (n=27) were retrieved from the Multi Organ Transplant Program at the Ajmera Transplant Centre (UHN). Serum samples were processed for targeted metabolomics processing by The Metabolomics Innovation Centre (TMIC Edmonton, AB, Canada, https://metabolomicscentre.ca) using a combination of direct injection mass spectrometry with a reverse-phase LC-MS/MS custom assay. This specialized assay enabled the accurate identification and quantification of 143 native metabolites encompassing amino acids, acyl carnitines, biogenic amines and their derivatives, uremic toxins, glycerophospholipids, sphingolipids, as well as various sugars. Mass spectrometric analysis was performed on an ABSciex 4000 Qtrap® tandem mass spectrometry instrument (Applied Biosystems/MDS Analytical Technologies, Foster City, CA) equipped with an Agilent 1260 or Waters series UHPLC system. The samples were delivered to the mass spectrometer by a LC method followed by a direct injection (DI) method.

2.2. Data Analysis

Metabolite concentrations reported in uM units (absolute concentrations) for each sample were used as input for MetaboAnalyst 5.0 software[22] for bioinformatics analysis. For each two-group comparison, the data underwent a series of processing steps, including: (i) removing features with more than 50% data missing, (ii) missing value imputation by replacing missing values with 1/5 of the minimum positive value, and (iii) normalization using quantile normalization, log2 transformation, and autoscaling. Subsequently, we applied multivariate analysis with Partial Least Square – Discriminant Analysis (PLS-DA) approach to identify significant metabolites based on Variable Importance in Projection (VIP) score calculated for each component. A metabolite was considered significant if VIP > 1. These metabolites were then categorized into their respective compound classes as per the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Human Metabolome Database (HMDB) using Pathway Analysis module of the software. Finally, we individually mapped them onto their primary biochemical pathways to gain a comprehensive visual representation of the metabolic changes. Boxplots illustrating the normalized expression of the significant metabolites and associated p-values from unpaired t-tests were generated using GraphPad Prism V.10.1.0 (GraphPad Software, San Diego, California, USA).

2.3. Integration of Clinical Variables with VIP Metabolites for Prediction of Patient Outcomes Post Transplantation

In this study, the cohort of 55 patients was divided into a train-test strategy with a split ratio of 75-25. Subsequently, on the training dataset, a feature selection step was performed to identify crucial metabolites exhibiting differential abundances across the three classes. This process involved receiver operating characteristic (ROC) curve analysis, performed using the filterVarImp function from the R package 'caret'[23]. Each metabolite underwent univariate evaluation, and three pairwise comparisons were conducted (Biliary vs MASH, Biliary vs TCMR, and MASH vs TCMR). The maximum area under the curve (AUC) was recorded for each pairwise comparison. Only metabolites surpassing an AUC threshold > 0.75 in at least two pairwise comparisons were considered significant. This feature selection procedure identified a subset of 20 metabolites deemed crucial for the classification task.

Subsequently, we integrated the metabolomic profiles with clinical and laboratory measurements of individuals to train a 3-way Random Forest classifier using the randomForest package[24] in R. Before model training, normalization procedures were applied to both the metabolomic profiles and clinical variables. To address class imbalances, particularly in the minority MASH and TCMR classes, and achieve a balanced class distribution, the SMOTE (Synthetic Minority Over-Sampling Technique) method[25] was employed to generate synthetic samples. The quality of samples generated using SMOTE was visually assessed through t-distributed stochastic neighbor embedding (t-SNE)[26] projections.

The efficacy of the three-class Random Forest classifier was assessed through the evaluation of the Out-of-Bag (OOB) error estimate. To further elucidate the multi-class classification problem, we employed a One-vs-Rest strategy, breaking it down into three distinct binary problems: Biliary versus Rest, MASH versus Rest, and TCMR versus Rest. Mitigating class imbalances in the binary outcomes involved under sampling the majority class. Moreover, we computed Area Under the Curve (AUC) values to gauge the effectiveness of each model in distinguishing between classes, leveraging a combination of metabolites and clinical markers.

To assess the significance of each explanatory variable (both metabolite and clinical variable), we employed a permutation-based variable-importance measure and Gini impurity criterion, which evaluated the capacity of predictors to mitigate data impurity or disorder. This assessment was conducted using out-of-bag data from the RF models.

3. Results

A total of 55 participants were enrolled in the study of which 10 had a diagnosis of post LT MASH, 18 had TCMR and 27 had biliary obstruction (Table 1). Most study participants were transplanted for steatotic liver disease and there was no significant demographic inter-group variability. Expectedly, alanine aminotransferase (ALT) levels were higher in participants with TCMR and MASH whilst alkaline phosphatase (ALP) levels were more elevated in individuals with biliary obstruction.

Within each two-group comparison, a total of 132 metabolites successfully passed the MetaboAnalyst processing steps detailed in the Materials and Methods section. Following these steps, the normalization was applied for further multivariate analysis.

3.1. Alanine, Aspartate, and Glutamate Metabolism Pathway Exhibited Notable Alterations in a Comparative Analysis of MASH (n=10) and TCMR (n=18) Patients

Employing PLS-DA, we identified 40 important metabolites, as documented in Table S1. The PLS-DA plot, as well as the top 15 features, are graphically depicted in Figure 1A, upper panel. Among these significant metabolites, several amino acids (serine, phenylalanine, alpha aminoadipic acid) and cholines (lysophosphatidylcholine acyl C18:1, lysophosphatidylcholine acyl C26:1) exhibited higher abundance in TCMR patients. Each of these amino acids have either been implicated in immune modulation or as biomarkers of steatohepatitis, in some cases both[27,28,29,30]. Serine deficiency, for example, has been repeatedly identified as a biomarker of MASH in non-transplant populations and seems to also have a well-defined role in T-cell responses[31].

The carnitines (nonaylcarnitine, decanoylcarnitine, octanoylcarnitine, dodecanoylcarnitine) were significantly more prevalent in the MASH group. Selected metabolites with their normalized concentrations are illustrated in Figure 1A, lower panel.

Conducting pathway analysis with VIP metabolites unveiled three significantly perturbed signaling pathways. Among these, the alanine, aspartate, and glutamate metabolism pathway was the most significantly affected, with an FDR-adjusted p-value of 1.09E-03 (Table 2).

3.2. Aminoacyl-tRNA Biosynthesis Emerged as the Pathway Exhibiting the Most Pronounced Alterations when Comparing MASH with Biliary Complications Group

As a secondary analysis, we undertook a comparison between MASH patients (n=10) and those with biliary complications (n=27). The application of Partial Least Squares Discriminant Analysis provided valuable insights, particularly through the visualization in the PLS-DA plot, which demonstrated a clear and effective separation of samples based on the respective diseases. To offer a more comprehensive understanding of the results, we have presented the top 15 metabolites in Figure 1B, upper panel.

From this analysis, we identified 39 metabolites with VIP scores greater than or equal to 1, and these findings are detailed in Table S2. Notably, our results indicated hydroxysphingomyeline as the top changed metabolite (Figure 1B, lower panel), with VIP > 2.9, most prevalent in MASH patients. Also, several phosphatidylcholines (lysophosphatidyl acyl cholines C14:0 and C18:0) and phosphatidylcholine acyl-alkyl C40:6 were more abundant in the MASH group. On the contrary, a two phosphatidyl diacyl cholines, C36:0 and C40:2, were found to be more prevalent in patients with biliary complications. Additionally, we noted an elevated abundance of serine and phenylalanine within this group. This differential metabolite distribution underscores the distinctive metabolic profiles associated with these distinct disease conditions. While some of these metabolites have not previously been identified in the development of either complication, a theoretical mechanism exists to explain their abundance or scarcity. For example, hydroxysphingomyeline C22: is a ceramide. Ceramide metabolism has been previously identified as being highly upregulated in non-transplant patients with MASH, supporting its potential as an effective biomarker in transplant populations[32].

Using the top 39 metabolites we identified three most significant metabolic pathways (Table 2), of which Aminoacyl-tRNA biosynthesis is the most significantly altered (FDR < 1.35E-06), with the contribution of eight metabolites, mostly amino acids (Asparagine, Phenylalanine, Glycine, Aspartate, Valine, Alanine, Isoleucine, Leucine, Tryptophan) from our list.

3.3. Significant Alterations in Butanoate (butyrate) Metabolism were Revealed in the Biliary vs TCMR Group Comparison

The third and final comparative analysis was performed between individuals with biliary complications and those in the TCMR group. A shorter list of just 28 significant metabolites (Table S3) was revealed by PLS-DA, the top 15 being illustrated in Figure 1C, upper panel. Within this set of metabolites, it was evident that serotonin and a pair of lysophosphatidyl acyl cholines (specifically C24:0 and C26:1) exhibited greater abundance in the TCMR group. In contrast, a couple of carnitines, specifically acetyl carnitine and pimeloyl carnitine, were more prevalent among patients with biliary complications (Figure 1C, lower panel). While serotonin has not previously been identified as a biomarker of post-LT complications, there is a strong biological mechanism for its marked difference in biliary vs TCMR groups, as serotonin is a potent modulator of T cells and has multiple functions in the liver including regulation of the biliary tree and cholangiocytes[33,34,35].

Furthermore, employing pathway analysis on metabolites with VIP scores exceeding 1, we identified three significantly altered metabolic pathways, as presented in Table 2. Among these pathways, the most pronounced alterations were observed in the butanoate (butyrate) metabolism pathway, characterized by a false discovery rate value of less than 8.13E-04.

In summary, the metabolomics analysis alone unveiled disease-specific metabolic alterations and highlighted the unique biological relevance of certain metabolites within these distinct pathways for each disease condition.

3.4. Integration of Clinical Information with Metabolomics Data

Prior to integration with conventional clinical and laboratory measurements, we performed feature selection on the metabolomics data using ROC curve analysis to retain only the most discriminative metabolites. Employing a random train-test split of 75%-25% samples, metabolite concentrations from the training set (n=40) were utilized to compute AUCs for each pairwise comparison, resulting in the elimination of 112 metabolites. Twenty metabolites were identified as crucial predictors, with AUC thresholds exceeding 0.75 for at least two of the classes. Notably, decanoylcarnitine C10, citric acid, succinic acid, phenylalanine, and serine emerged as the top-ranking metabolites, and their details are presented in Table 3.

Subsequently, we developed a Random Forest classifier model by combining the eight key clinical variables: age, sex, alkaline phosphatase (ALP), alanine transaminase (ALT), aspartate aminotransferase (AST), creatinine, hemoglobin (HGB), and primary indication for transplantation—integrated with the top 20 selected metabolites. The t-SNE projections of the original dataset are shown in Figure S1A. The SMOTE-generated samples have good coherence with the original samples as shown in Figure S1B. A 3-way classification was performed to predict the likelihood of a patient belonging to one of three classes: Biliary, MASH, or TCMR. Model parameters were tuned using the out-of-bag (OOB) error, revealing 5 as the optimal number of candidate predictors randomly drawn for a split and the number of trees as 500, based on minimum OOB error. The 3-class classification model yielded an overall OOB estimate of the error rate at 19.75%.

Figure 2A illustrates that at the individual class level, the model demonstrated the maximum ability to distinguish MASH samples with an OOB error estimate of only 7.4%, compared to 22.2% and 29.6% for Biliary and TCMR groups, respectively. The 3-way classifier model achieves an overall accuracy of 79.66%. Serotonin and serine have surfaced as the primary predictors, identified through both Mean Decrease in Accuracy and assessment via the Gini Impurity Criterion, which gauges the capacity of predictors to mitigate data impurity or disorder.

Apart from the 3-way classification, we also evaluated three Random Forest models predicting binary outcomes: Biliary model (Biliary vs. Rest), MASH model (MASH vs. Rest), and TCMR model (TCMR vs. Rest). The resulting OOB error rates were found to be 22.64%, 5.66%, and 24.53% for the Biliary, MASH, and TCMR models, respectively. Figure 2 (B—D) show the ROC curves and the rank of the variables for the three models. The AUCs for the Biliary, MASH, and TCMR models are 0.882, 0.972, and 0.896, respectively.

In our variable importance analysis, as shown in Figure 3, serotonin emerged as a top predictor for the 3-way-Classification, Biliary, and TCMR models, consistent with Mean Decrease in Accuracy and Mean Decrease in GINI.

The amino acid serine and hydroxysphingomyeline (SM(OH)C22:1) emerged as leading predictors for 3-way-Classification and the MASH model, indicating their importance in classifying MASH samples. Additionally, phenylalanine, decanoylcarnitine, and kynurenine were crucial predictors of MASH. The liver enzymes AST and ALT appeared as the topmost important clinical variables. Abundance levels of top predictor metabolites are plotted in Figure 4.

3.5. Comparative Analysis of Integrated Model versus Individual Modalities

We additionally assessed our integrated 3-way classification model, which combines clinical variables and metabolites, alongside two other 3-way classification models trained solely on a single data type: i) clinical variables alone, and ii) metabolites alone. Our Random Forest classifier, when trained on solely clinical variables, produced an overall OOB estimate of the error rate at 25.93%. At the individual class level, depicted in Figure 5A, the model distinguished both MASH and TCMR groups, with OOB error estimates of 22.2% and 33.3% for the Biliary group. The 3-way clinical-only model achieves an overall accuracy of 73.37%. The top clinical predictors were ALT and AST as shown in Figure 5B. While these conventional clinical variables can indicate liver graft injury, they are not specific to etiologies of graft pathology.

In contrast, our Random Forest classifier trained on top ranked metabolites alone, produced a lower overall OOB estimate of the error rate at 22.22%. Depicted in Figure 5C, the confusion matrix for the metabolites-only model illustrated improved performance at individual class levels, particularly evident in the MASH group with an OOB error estimate of 14.81% compared to the clinical-only model. Achieving an overall accuracy of 77.14%, the metabolite-only model identified Serine and Serotonin as top predictor metabolites, based on the mean decrease in accuracy, as shown in Figure 5D.

4. Discussion

In this study, we aimed to identify distinct, measurable metabolomic profiles to differentiate causes of liver graft injury along with clinical variables in post-transplant populations. Our metabolomic analysis detected specific metabolites that exhibited significant changes in individuals with TCMR, biliary complications and MASH post-liver transplant. In addition to identifying individual metabolites that varied between disease groups, pathway analysis was conducted to determine which metabolic pathways were most differentially affected.

In the case of MASH vs TCMR, amino acid metabolism was broadly highlighted, with two of the main pathways affected being alanine aspartate glutamate metabolism and arginine biosynthesis. Both pathways have been affected in steatohepatitis in non-transplant populations and in rejection for transplant populations, further supporting their potential inclusion as relevant biomarkers in post-LT complications[36,37,38].

Comparing MASH vs Biliary groups, amino acid metabolism was once again of great significance, with branch chain AA(BCAA) synthesis and alanine aspartate glutamate metabolism being highlighted. The BCAA synthesis pathway may be particularly useful as a biomarker of MASH, as levels of these amino acids have been shown to increase greatly in non-transplant MASH populations[39].

Lastly, the biliary vs TCMR groups highlighted one unique metabolic pathway, butanoate metabolism. Butanoate is a short chain fatty acid synthesized by enteric bacteria and has been repeatedly implicated in mediating tolerogenic phenotypes of T cells and may therefore be associated with graft tolerance or rejection[40,41]. Considering its unique role in mediating immune functioning, biomarkers associated with butanoate metabolism may be effective in identifying TCMR.

We then developed an ML tool, a Random Forest classifier to predict liver graft pathology and compared its accuracy on three circumstances: when using only (i) clinical parameters, (ii) metabolites, and (iii) integrating both categories. The improvement in overall accuracy observed in the integrated model, as discussed in Section 3.4, reaffirms our assertion that the complete complexity of graft pathology cannot be adequately captured by any singular data modality. Our results show that clinical variables and metabolites are complementary in nature and integrating them provides a more comprehensive understanding of liver graft pathology.

When integrating metabolomic data with clinical information via ML modeling, serotonin was identified as a top predictor for the 3-way-Classification. Serotonin has been shown to regulate fibrosis progression[42,43], with higher levels being linked to increased chances of developing MASLD[44]. Our ML model also revealed serine and hydroxysphingomyeline (SM(OH)C22:1) as crucial metabolites for 3-way-Classification and the MASH model. Indeed, previous studies have shown the association between serine deficiency and MASH[45,46]. Sphingomyelin is the most frequently observed sphingolipid in mammalian cells and circulating levels have shown promise as a noninvasive biomarker of MASH[45]. It is considered a bioactive lipid, functioning as a component of cell membranes, with a role in cell signaling, growth, death, senescence, adhesion and migration. The liver plays an important role in lipid metabolism, taking up free fatty acids, forming triglycerides and VLDL in addition to acting as a site for lipid storage. Lovric et al identified a positive correlation between serum concentrations of sphingomyelin with increased ectopic fat accumulation (including hepatic steatosis)[47,48]. In murine models it has been shown that sphingomyeline synthesis is activated in MASH, being related to hepatocyte pyroptosis[49]. However, Zhou et al noted a significant decrease in serum sphingomyelin and lysophosphatidylcholine levels in individuals with MASH versus MASLD[50].

Other top metabolites as predictors for MASH were phenylalanine, decanoylcarnitine, and kynurenine. In agreement with our results, other studies have indicated an increased level of the phenylalanine in MASH patients[51]. Hanssen et al. have shown alteration of the kynurenine pathway in MASLD patients, favoring inflammation and fibrosis through regulation by inflammation markers such as IFN, IL6, LPS[52].

Interestingly, our results showed that carnitines were more prevalent in MASH. Studies looking at the different diets on the serum metabolomic profile noted an increase in both long and short chain acyl carnitines with a traditionally more obeso-genic, ‘western’ diet when compared to individuals on a vegetarian diet and associated with increased risk of cardiovascular disease[42,43]. Additionally, there was a positive correlation between short chain acyl carnitines and fasting insulin levels[42]. Under normal physiologic conditions, carnitine plays an important role in the ß oxidation of fatty acids by facilitating the transmembrane transfer of acetyl-CoA via acetyl carnitine. In pathological situations, including obesity and MASH, where there is an excess of fatty acids due to insulin resistance, increased lipogenesis and impaired fatty acid oxidation, acyl carnitine levels have been shown to accumulate[45,46], with one study demonstrating an AUC > 0.90 for MASH with elevated levels of long chain acyl carnitines: C20, C16:1 and C14:1OH, in the pre-transplant population[45].

5. Conclusions

In this project, we have developed an ML tool integrating serum metabolites with clinical variables in liver transplant patients with MASH, TCMR and biliary complications. Our tool appears to be a promising non-invasive indicator for detecting graft pathology. The model identified serine and serotonin as top altered metabolites, and liver enzymes AST and ALT as the most important clinical variables. It also exceled in predicting the occurrence of MASH following a transplant with the highest accuracy, with an OOB error estimate of 7.4% compared to 22.2% for biliary and 29.6% for TCMR. In the prediction of binary outcomes across three models: Biliary (biliary vs. others), MASH (MASH vs. others), and TCMR (TCMR vs. others), the Area Under the Curve (AUC) scores were 0.882, 0.972, and 0.96, respectively.

As a limitation to note, our dataset accurately reflects the group of patients from our own institution, without incorporating data from an external cohort for validation purposes. Despite this limitation, it is important to recognize that this work serves as a pilot study featuring a thoroughly detailed and carefully selected group of patients. The outcomes from this study are promising, highlighting the potential of our approach.