Serum metabolite profiles are associated with the presence of advanced liver fibrosis in Chinese patients with chronic hepatitis B viral infection

Study design and participants

Two data sets were enrolled in this study. Cohort 1 was recruited between April 2013 and June 2015 at Shuguang Hospital Affiliated to Shanghai University of Traditional Chinese Medicine, consisted of 1,006 participants, including 504 CLD patients with chronic HBV infection and 502 NC as our training cohort to identify serum metabolite markers and establish predictive models (Table 1). More detailed inclusion and exclusion criteria can be found in the Supporting Information.

Cohort 2, recruited between December 2016 and December 2017 at Xiamen Hospital of Traditional Chinese Medicine, consisted of 300 CLD patients with chronic HBV infection and 90 NC. Data obtained from Cohort 2 were used as a validation set to further verify the performance of the models established from the Cohort 1. Detailed information about this cohort can be found in the Supporting Information. Sample size was not determined by statistical methods and was comparable to other studies in the field [4-8, 21, 22].

In this study, the diagnosis and the sample collection were performed using exactly the same protocols to avoid “external” influences. The samples were provided to lab staffs blind samples with respect to patient identity and other clinical information.

The study was approved by the institutional review board of each hospital. All participants provided written informed consent.

Classification Performance Evaluation

ROC curve is a plot of the true positive rate (sensitivity/recall) against the false positive rate (1-specificity) at different cutoffs of a binary classifier. AUROC measures the area under the ROC curves and a higher value of AUROC suggests better classification performances while an AUROC of 0.5 represents the random guess. The PR curve demonstrates the relationship between positive predictive values (precision) and true positive rate (sensitivity/recall), and a higher value of AUPR indicates better diagnostic capacity of the model. PR curves are usually preferable for evaluating unbalanced data compared to ROC curves. NRI and IDI were also used for the evaluation of prediction improvement. We compared RF models to existing clinical indices by splitting the continuous risk scores into ten equal risk intervals (default). We used the R software version 3.2.3 for data analysis and the “PRROC” R package for binary ROC and PR curves[23], the “pROC” package for calculating the specificities and sensitivities of classifiers[24], and the “PredictABEL” package for NRI and IDI calculation[25].

Feature Selection and Method Comparison

Quantitative variables were expressed as mean ± SD for clinical parameters and median (25% quantile, 75% quantile) of log10 transformed concentration for metabolites. Categorical variables were expressed as percentages. The univariate analysis (Wilcoxon rank-sum test) was carried out to identify the variables that were significantly different between CLD patients and NC, between fibrosis and cirrhosis (S0-3 vs. S4), and among CLD patients at different fibrotic stages (early stage fibrosis (S0-2) vs. advanced stage fibrosis (S3-4)).

For differential metabolites with p < 0.001 across all univariate analyses were used in two machine learning methods, LASSO[26] and RF[27] to further select markers for the three classifications listed above. Data were log and z-score transformed before being fed into LASSO to ensure that the coefficients were comparable with each other. The regularization parameter lambda of LASSO was determined using 10-fold CV. The RF model used 500 decision trees. We ranked the metabolites according to their LASSO non-zero coefficients and RF mean decrease of accuracy, and kept the intersection of top 5 LASSO and RF metabolites in the three classifications. Considering the overlaps of the second and the third classification tasks, we further selected the intersecting variables of these two situations and then, merged with variables selected from the first situation to construct our final metabolite markers (Fig. 2).

• Download figure
• Open in new tab

Figure 2. Workflow chart of feature selection.

For a total of 98 metabolites (including AAs, BAs and FFAs), univariate analyses (Wilcoxon rank-sum test) were employed for three clinical aims (aim-1: CLD vs. NC, aim-2: Cirrhosis vs. Fibrosis, aim-3: Early fibrosis vs. Advanced fibrosis). 26 metabolites with P < 0.001 in all three clinical aims were selected and fed into LASSO and Random Forest algorithms for three aims. The overlap of top-5 LASSO non-zero coefficients and top-5 important variables from Random Forest (by mean decrease of accuracy) were selected. For aim-2 and aim-3, we selected the overlapped variables and combined with variables selected from aim-1 to yield the final panel four metabolites.

Note: “OR” means the union of two sets, “AND” means the intersection of two or more sets.

To identify an appropriate classification method, we introduced two linear models, i.e., logistic regression (LR), linear discriminant analysis (LDA) and one decision tree-based ensemble model, i.e., RF, for the classifier construction for the markers we selected. For RF, we used 500 decision trees and two candidate variables at each split. For LDA, the tolerance parameter was set to 1.0E-4 (default). We applied 10-fold CV on the training set (Cohort 1) to compare the classification performances of these four models and three established fibrosis markers, i.e., APRI, AST/ALT ratio, and FIB-4. AUROC and AUPR were recorded at each internal validation set in CV. We used R packages “randomForest”, “glmnet” and “MASS” for RF, LASSO and LDA constructions respectively [28, 29].

Predictive Model Construction and Validation

We trained the final RF models for different classification objectives using the training data (Cohort 1), with Model 1 differentiating CLD and NC, Model 2 differentiating fibrosis and cirrhosis, and Model 3 differentiating early and advanced stages of liver fibrosis. A total of 500 decision trees were included in a single RF model with two variables randomly sampled as candidates at each split. We re-balanced the sample size for different groups at each bootstrap resampling step for Models 2 and 3 considering the unbalanced samples.[30]

In RF, each decision tree was fitted on the bootstrap samples and tested on the untouched OOB-samples. Thus, the OOB-predictions provided unbiased estimates of how the RF model performed on the training data and were used for the evaluation on Cohort 1. We further validated our mark panel-based RF models in the independent validation data sets from Cohort 2, and compared results with the established fibrosis markers, AST/ALT ratio, APRI and FIB-4. ROC and PR curves were drawn and AUROC and AUPR values, respectively, were calculated to evaluate their diagnostic performances. Optimal cut-offs were selected to maximize the sum of sensitivity and specificity for the RF model. For APRI, FIB-4 and AST/ALT, predefined cut-offs were used (1.0 and 2.0 for APRI to distinguish fibrosis and cirrhosis,[6] 1.45 and 3.25 for FIB-4 to distinguish S0-2 and S3-4,[7] and 0.8 and 1.0 for AST/ALT to distinguish S0-2 and S3-4[5, 31]). Bootstrap resampling (1,000 times) was conducted to calculate 95% confidence intervals (CIs) of AUCs for all binary classifiers. A comparison of the AUROC of our biomarker panel vs. FIB-4, AST/ALT, or APRI was performed using DeLong’s test. The significance level was adjusted for multiple testing according to the Benjamin and Hochberg procedure.[32] Log and z-score transformed data were also used for constructing heatmaps. The R packages “ggplot2” and “cowplot” were used for data visualization and multiple plots arrangement.

We further derived an RF risk score for each participant based on the marker panel and logit function of the predicted probability (Prob.) of the RF model for corresponding classification objective:

F1 scores were then calculated at the predefined cutoffs using following formula:

To determine whether the RF score could independently predict the fibrosis staging in the presence of other potential confounding factors, we applied logistic regression on the RF-score to differentiate cirrhosis from fibrosis as well as discriminate early and advance fibrosis while adjusting for HBV DNA levels, the degree of necro-inflammation, HBeAb status, HBeAg status, liver function tests (i.e., PT, ALB, DBIL, IBIL), platelets, BMI, and medication (entecavir) use.

Multi-group Classification of S0-2 vs. S3 vs. S4

We built a new RF model based on our metabolite marker panel and applied multinomial regressions to APRI, AST/ALT, and FIB-4 separately to differentiate S0-2 vs. S3 vs. S4 in Cohort 1. Then, we compared and validated these multi-group classifiers on both Cohort 1 and Cohort 2 data sets using micro-average ROC and PR curves. Micro-average ROC and PR curves were calculated by stacking binary classification results from each group together to generate a concatenated binary classification result.[33] We then calculated AUROC and AUPR with 95% CIs using 100 times bootstrap resampling. We used the “multiROC” R package for calculating the micro-average AUROC and AUPR as well as for plotting[34].

Code availability

R, GraphPad, and in-house metabolomic processing pipeline were used for the data analyses presented in this study. We used GraphPad Prism 7, R software version 3.2.3 with packages including PRROC, pROC, PredictABEL, randomForest, glmnet, MASS, ggplot2, cowplot and multiROC. The metabolomics data were processed and analyzed by TargetLynx (Waters). The computer code used to generate the results reported in this study is available from the authors upon request.

Data availability statement

All the data supporting the findings of this study are available within the article and its Supplementary Information files or from the corresponding author upon reasonable request.

Characteristics of the participants

Two independent cohorts were studied (Fig. 1). CLD groups were staged and assigned according to the results of their liver biopsy. Cohort 1 consisted of 1,006 participants (502 NC and 504 biopsy-proven HBV-CLD patients (400 with liver fibrosis (S0-3) and 104 with cirrhosis (S4); or 349 with early stage fibrosis (S0-2) and 155 with advanced stage fibrosis (S3-4)). Cohort 2 consisted of 390 participants (90 NC and 300 biopsy-proven CLD patients comprising 141 with fibrosis and 159 with cirrhosis, or 134 with early stage fibrosis (S0-2) and 166 with advanced stage fibrosis (S3-4). Models established from Cohort 1 were validated in Cohort 2 (Models 1, 2 and 3). The cohort stratification and major demographic and clinical characteristics are shown in Table 1. More detailed clinical data are provided in Table S1.

• Download figure
• Open in new tab

Figure 1. Study design.

Serum metabolites were quantified in Cohort 1 (504 biopsy-proven HBV-CLD patients and 502 NC) and were used to identify candidate markers. After data analysis and feature selection, four metabolites were selected to compose our marker panel. Different machine models and clinical indices were compared using 10-fold cross validation. Three RF models were constructed to diagnose CLD from NC (Model 1), differentiate fibrosis vs. cirrhosis (Model 2) and grade early fibrosis vs. advanced fibrosis (Model 3) in Cohort 1. These three were further validated in the independent HBV Cohort 2.

Quantification of metabolites in serum

Using targeted metabolomic protocols establlished in our lab,[35-37] we quantified the concentrations of 98 metabolites, including 24 BAs, 42 FFAs and 32 AAs, in the sera of all participants (Table S2). These metabolites were used for the subsequent metabolite marker selection.

Serum metabolite marker selection

From the 98 serum metabolites, we identified 26 differential metabolites in three classification situations (i.e., to diagnose CLD patients from NC, to differentiate fibrosis from cirrhosis, and to differentiate advanced fibrosis from early fibrosis) using univariate analysis (Wilcoxon rank sum test, P<0.001). The 26 statistically significant metabolites were then entered into least absolute shrinkage and selection operator (LASSO)[26] and random forest (RF)[27]. According to the rank of LASSO non-zero coefficients and RF mean decrease of accuracy, four metabolite markers were selected, which included one FFA, linoelaidic acid (C18:2 n6t), one BA, taurocholate (TCA), and two AAs, tyrosine (Tyr) and valine (Val) (Fig. 2). The principal component analysis (PCA) of these four metabolite markers also showed a clear separation between CLD patients and NC (Fig. S1). We also derived one ratio, the Tyr/Val ratio to further improve the classification performances while also including one extra accessible risk factor, age, to our panel for the differentiation of fibrosis and cirrhosis and the staging of fibrosis. Correlations of the four metabolites with fibrosis stage, necro-inflammation, CPA, AST, ALT, AST/ALT ratio, PLT, FIB-4, and APRI were assessed using Spearman correlation analysis (Fig. S2). The four metabolite markers (including the Tyr/Val ratio) all significantly correlated with fibrosis stage (ρ=0.38 for TCA, ρ=0.50 for Tyr, ρ=0.53 for Tyr/Val ratio and ρ=0.23 for C18:2 n6t) using Spearman’s correlation analysis. In addition, we found our metabolite markers showed stronger associations with the fibrosis stage than the previously used clinical indices.

To determine an appropriate classification model, we applied 10-fold cross-validation (CV) on Cohort 1 to compare the classification performances of RF models and two linear models (i.e., logistic regression (LR), linear discriminant analysis (LDA)) as well as the clinical indices, APRI, AST/ALT ratio, and FIB-4. The CV-area under the receiver operating characteristic curve (CV-AUROC) and area under the precision-recall curve (CV-AUPR) were employed as the evaluation metrics. We found that, to differentiate CLD from control, APRI, LR, LDA and RF had the highest AUROCs and AUPR, while RF demonstrated the most robust classification performance (Fig. S3a). For the differentiation of fibrosis and cirrhosis and S0-2 vs. S3-4, RF outperformed other methods with the highest CV-AUROC and CV-AUPR overall (Figs. S3b, c). PCA scores plot showed linearly separable discrimination between the most CLD and control subjects (Fig. S1), thus linear models could achieve good classification performances. However, for a situation where there is more extensive overlapping of groups (Fig. S4), the decision tree-based ensemble learning algorithm, RF, achieved improved classification performances compared to other methods (Fig. S3).

Model 1: Differentiating CLD patients from NC

The concentration of linoelaidic acid (C18:2 n6t) was significantly higher in the control group than in the CLD group, conversely the levels of TCA, and Tyr, and Tyr/Val ratio were higher in the CLD group than in the control group (Figs. 3a, b).

• Download figure
• Open in new tab

Figure 3. Metabolite marker panel and Model 1 for CLD with chronic HBV infection diagnosis.

(a) Comparison of the four markers between CLD patients and NC in Cohorts 1 and 2. (b) Waterfall plot of RF-score and corresponding heatmap for the four markers in all data sets. (c) ROC curves of Model 1 (RF model constructed with four markers), APRI, AST/ALT and FIB-4 in Cohort 1. (d) PR curves of Model 1, APRI, AST/ALT and FIB-4 in Cohort 1. (e) ROC curves of Model 1, APRI, AST/ALT and FIB-4 in Cohort 2 validation set. (f) PR curves of Model 1, APRI, AST/ALT and FIB-4 in Cohort 2 validation set. (g) The diagnosis RF-score in NC and CLD patients in training and validation sets. *** p < 0.001, Wilcoxon rank sum test. The optimal cut-off value of the RF-score was 0.434.

Model 1 was constructed using an RF model that utilized these four metabolite markers, to differentiate CLD patients from NC in Cohort 1. Out-of-bag (OOB) estimates were employed for the RF model evaluations. Model 1 showed an AUROC of 0.997 (0.993-1.000) and AUPR of 0.994 (0.998-1.000) (Figs. 3c, d) which was significantly higher than the APRI (AUROC=0.973, p<0.001), FIB-4 (AUROC=0.848, p<0.001) and AST/ALT ratio (AUROC=0.665, p<0.001) (Table 2). An example decision tree from the RF model is shown in Fig. S5a, where we observed that the lower concentration of C18:2 n6t and the higher concentration of TCA would lead to higher risk of CLD.

Based on the OOB predicted probabilities, we calculated a diagnostic RF-score for Model 1 using the logit function. The waterfall plot showed a clear ascending trend of RF-scores from NC (lower RF-scores) to CLD patients (higher RF scores) along with the differentiation trend shown in the heatmap of the four markers (Fig. 3b). We observed significant differences in the RF-score between both groups in Cohort 1 (p<0.001, Fig. 3g), yielding a sensitivity of 98.4% and specificity of 99% for CLD patients in the training set at a cutoff value of 0.434 (Table 2). The sensitivity and specificity of our RF model were superior to those of AST/ALT ratio, APRI, and FIB-4 for differentiating CLD patients from NC using the optimal cutoffs generated in Cohort 1 using the Youden index (Table S4).

Model 2: Differentiating cirrhosis from fibrosis among CLD patients

The discriminant prediction model was constructed using an RF model employing the four metabolite markers along with age to differentiate CLD patients with cirrhosis from those without cirrhosis in Cohort 1. This model demonstrated an AUROC of 0.941 (0.914-0.964) and AUPR of 0.87 (0.824-0.913) (Figs. 4a, b) based on OOB predictions. These results were better than those of the APRI (AUROC=0.698, p<0.001), AST/ALT (AUROC=0.815, p<0.001) and FIB-4 (AUROC=0.869, p<0.001) (Table 2). We showed an example decision tree for Model 2 in Fig. S5b, and we found that the higher the Tyr/Val ratio, Tyr and C18:2 n6t, the higher the risk of CLD with cirrhosis.

• Download figure
• Open in new tab

Figure 4. Model 2 for differentiating fibrosis vs. cirrhosis and Model 3 for differentiating early fibrosis vs. advanced fibrosis in CLD patients with chronic HBV infection.

(a) ROC curves of Model 2 (RF model constructed with four metabolite markers and age), APRI, AST/ALT, and FIB-4 in Cohort 1. (b) PR curves of Model 2, APRI, AST/ALT, and FIB-4 in Cohort 1. (c) ROC curves for Model 2, APRI, AST/ALT, and FIB-4 for the Cohort 2 validation set. (d) PR curves for Model 2, APRI, AST/ALT, and FIB-4 for the Cohort 2 validation set. (e) The RF-score in CLD patients with fibrosis and cirrhosis in the HBV training, validation sets. The optimal cut-off value of the RF-score was 0.01. (f) ROC curves of Model 3 (RF model constructed with four metabolite markers and age), APRI, AST/ALT, and FIB-4 in Cohort 1. (b) PR curves of Model 3, APRI, AST/ALT, and FIB-4 in Cohort 1. (c) ROC curves for Model 3, APRI, AST/ALT, and FIB-4 for the Cohort 2 validation set. (d) PR curves for Model 3, APRI, AST/ALT, and FIB-4 for the Cohort 2 validation set. (e) The RF-score in CLD patients with S0-2 and S3-4 in the HBV training, validation sets. The optimal cut-off value of the RF-score was −0.115. *** p < 0.001, Wilcoxon rank sum test.

The Model 2 RF-score differentiated CLD patients with cirrhosis from fibrosis in Cohort 1 (p<0.001) (Fig. 4e). The constructed model yielded a sensitivity of 87.0% and specificity of 90.4% in the Cohort 1 data set at a cutoff value of 0.01 (Table 2). The RF-scores remained significant with a coefficient of 0.755 (p<0.001) after adjusting for HBV DNA levels, degree of necro-inflammation, HBeAb status, HBeAg status, body mass index (BMI), platelets (PLT), liver function tests (i.e., prothrombin time (PT), albumin (ALB), direct bilirubin (DBIL), indirect bilirubin (IBIL)), and medication (Entecavir) (Table S3). The accuracy of our RF model was superior to those of AST/ALT ratio, APRI, and FIB-4 (Table 2 and Table S4).

Model 3: Differentiating advanced fibrosis from early fibrosis among CLD patients

In this study, fibrosis stages 0-2 were defined as early fibrosis, and stages 3-4 were defined as advanced fibrosis. Model 3 was established based on age and the four metabolite markers selected from Cohort 1 data using the RF model. It was then shown to successfully separate CLD patients with early fibrosis from those with advanced fibrosis in Cohort 1 with AUROC of 0.918 (0.889-0.946) and AUPR = 0.892 (0.854 - 0.925) (Figs. 4f, g). Model 3 results demonstrated better classification performances than those of APRI (AUROC=0.647, p<0.001), AST/ALT (AUROC=0.714, p<0.001) and FIB-4 (AUROC=0.802, p<0.001) in predicting liver fibrosis stages (Table 2). An example decision tree from Model 3 showed that the higher Tyr/Val ratio, Tyr, age and TCA indicated a higher risk of CLD with advanced firbosis (Fig. S5c).

A logit diagnostic RF-score for Model 3 differentiated CLD patients with early stage fibrosis from those with advanced fibrosis in Cohort 1 (Fig. 4j). The model yielded a sensitivity of 86.7% and specificity of 90.5% in Cohort 1 at a cutoff value of −0.115 (Table 2). After adjusting for HBV DNA levels, degree of necro-inflammation, HBeAb status, HBeAg status, liver function tests (i.e., PT, ALB, DBIL, IBIL), platelets, BMI, and medication (Entecavir) use, RF-scores remained statistically significant with a coefficient of 0.805 (p<0.001) (Table S3). The accuracy of our RF model was superior to those of AST/ALT ratio, APRI, and FIB-4 (Table 2 and Table S4).

Validation of the predictive models in an independent HBV cohort (Cohort 2)

The metabolite markers identified and related models obtained in Cohort 1 were further validated for their liver fibrosis staging performance as well as for CLD diagnosis performance in Cohort 2, and the results were similar to those obtained from Cohort 1 (Table 2).

For the diagnosis of CLD patients, compared to APRI (AUROC = 0.879, AUPR = 0.958), AST/ALT (AUROC = 0.603, AUPR = 0.849), and FIB-4 (AUROC = 0.707, AUPR = 0.897), we again observed higher classification performances for Model 1 with AUROC of 0.977 (0.963-0.988) and AUPR of 0.993 (0.989-0.997) in the validation set (Figs. 3e, f). In addition, the Model 1 predicted RF-score in Cohort 2 differentiated CLD from NC with a sensitivity of 92.2% and specificity of 94.4% at the cutoff values determined for Cohort 1 (Fig. 3g, Table 2).

Applying Model 2 to Cohort 2 successfully discriminated cirrhotic patients from fibrotic patients with an AUROC of 0.844 (0.797-0.884) and AUPR of 0.827 (0.761-0.884) (Figs. 4c, d) and outperformed those of the APRI (AUROC=0.608, p<0.001), AST/ALT (AUROC=0.684, p<0.001) and FIB-4 (AUROC=0.758, p<0.001) indices. The Model 2 RF-score in Cohort 2 differentiated cirrhotic patients from fibrotic patients with a sensitivity of 71.8% and specificity of 81.6% at the same cutoff value used for the Cohort 1 data set (Fig. 4e, Table 2). Similarly, applying Model 3 to grade fibrosis stage in Cohort 2 resulted in greater performance with AUROC of 0.807 (0.756-0.852) and AUPR of 0.817 (0.764-0.866) (Figs. 4h, i) than those of APRI (AUROC=0.595, p<0.001), AST/ALT (AUROC=0.667, p<0.001) and FIB-4 (AUROC=0.739, p=0.01) indices. And the Model 3 RF-score also differentiated S0-2 fibrosis from S3-4 fibrosis with a sensitivity of 72.9% and specificity of 76.1% (Fig. 4j, Table 2).

We then introduced net reclassification improvement (NRI) and integral discriminant improvement (IDI) to quantify the improvement of our model to existing clinical indices. For different classification aims (Control vs. CLD, Fibrosis vs. Cirrhosis, S0-2 vs. S3-4) in an independent validation cohort (Cohort 2), the categorical and the continuous NRI and IDI of the RF models all achieved positive values when compared to FIB-4, APRI and AST/ALT, suggesting an augmentation of classification performances for our biomarker panel and RF models (Table S5).

Classification of S0-2 vs. S3 vs. S4

In addition to the binary classifications that we have performed, we further determined whether our biomarker panel could classify multiple groups among CLD patients. We trained a new RF model with the marker panel and applied multinomial regression to APRI, AST/ALT and FIB-4 respectively for the discrimination of S0-2 vs. S3 vs. S4 using Cohort 1. We compared their performances on Cohort 1 (OOB predictions of RF model) and Cohort 2 using micro-average AUROC and AUPR and we found that our marker panel-based multi-group classifier outperformed other methods. In the Cohort 1 data, our classifier showed higher AUROC of 0.944 (0.928-0.963) and AUPR of 0.908 (0.883-0.938) compared to APRI (AUROC=0.79, AUPR=0.658), AST/ALT (AUROC=0.817, AUPR=0.688), and FIB-4 (AUROC=0.858, AUPR=0.774) (Fig. S6a, b). In the Cohort 2 validation data, our marker panel classifier consistently displayed higher AUROC of 0.841 (0.799-0.885) and AUPR of 0.748 (0.674-0.81) compared to APRI (AUROC=0.790, AUPR=0.608), AST/ALT (AUROC=0.772, AUPR=0.597), and FIB-4 (AUROC=0.816, AUPR=0.699) (Figs. S6c, d).