ATB-DILI longitudinal cohort: a prospective study with high-density sampling

We enrolled 118 newly diagnosed TB patients and established a longitudinal cohort with multiple timepoints blood sampling (baseline and 1 W, 2 W, 3 W, 4 W, 6 W, 8 W after medication). A total of 14 individuals who developed DILI were included in ATB-DILI group and age-, gender-, and body mass index (BMI)-matched 14 individuals who complete their treatment without any adverse drug effects were included in ATB-Ctrl group.

Demographics and clinical baseline characteristics are presented in Table 1 and therapeutic regimes for each patient are listed in the supplementary Table 7. Dynamic changes of ALT and AST were displayed in the supplementary Fig. 1. ATB-DILI subjects were mainly presented with hepatocellular hepatitis that was characterized by ALT and AST elevation and the onset time of ATB-DILI was different due to individual differences. Based on sample collection time, 9 individuals of each group who were enrolled in the study from November 2019 to December 2020 were assigned to a discovery cohort and age-, gender-, and BMI-matched 18 healthy controls were also involved in the cohort. (Fig. 1A, supplementary Table 8).

Table 1 The baseline clinical features of the longitudinal cohortFig. 1

Untargeted metabolomics reveals dynamics changes of ATB-DILI. (A): The ATB-DILI group: T1 (baseline), T2 (one-third course), T3 (two-third course), T4 (ATB-DILI), and T5 (recovery period). The ATB-Ctrl group: T1 (baseline), T2 (2 W), T3 (4 W), T4 (6 W), and T5 (8 W). (B) Metabolites compositions of global metabolomics. (C) Sample distribution of PLS-DA. (D) Heatmap displayed the average intensity of top77 metabolites

ATB-DILI progression is associated with dynamic metabolic changes

A total of 107 samples within the discovery cohort were analyzed for global metabolomics. After quality control and peak identification, we successfully identified 895 metabolites across the different samples. Of these, lipid metabolites account for 46.82% (419/895) which is dominated by fatty acids (FAs) and bile acids (BAs) (48.93%, 205/419) (Fig. 1B). Then, we examined the data globally with PLS-DA, in which the samples were discriminated from the different three groups based on the first two components (Fig. 1C; the PCA results in supplementary Fig. 2).

Next, we performed a comparison of different groups to find statistically significant changed metabolites and data revealed the number of differential metabolites increased along with disease progression in the ATB-DILI group (Supplementary Figs. 3, 4, 5, 6, 7). To identify potential predictive biomarkers of ATB-DILI, we screened the differential metabolites using the following criteria: (1) ATB-DILI T4/T1: fold change (FC) ≥ 1.5 or FC ≤ 0.66, FDR ˂ 0.05, (2) ATB-DILI T4/ATB-Ctrl T4: FC ≥ 1.5 or FC ≤ 0.66, (3) exclusion of exogenous drugs. A total of 77 compounds were obtained, of which 62 metabolites significantly increased at T4 while 15 metabolites decreased. Hierarchical clustering based on the top 77 metabolites revealed that FAs significantly increased at T2, reached their peak at T4, and decreased at T5 while BAs elevated since T3, reached their peak at T4, and remained at a high level at T5 (Fig. 1D).

Metabolite groups altered during ATB-DILI

We further explored the expression characteristics of different subgroups metabolites. The long-chain saturated and unsaturated FAs exhibited a gradual increase at T2-T4 in the ATB-DILI group, with levels higher than those observed in the ATB-Ctrl and HC groups (Fig. 2A). Concurrently, long-chain carnitines that are responsible for FAs transportation from cytoplasm to mitochondrion had also increased at T2-T4 in the ATB-DILI group, and ketone body also exhibited a similar expression pattern (Fig. 2B and C). The phenomenon suggests that FAs transportation, β-oxidation, and ketone body synthesis are enhanced as the disease progresses. Furthermore, we also noted a significant increase in various BAs at T3-T4 in the ATB-DILI group (Fig. 2D), indicating a disturbance in the metabolic balance of BAs as the disease progresses.

Fig. 2

Functional metabolites alter during ATB-DILI. Expression characteristics of long-chain FAs (A), long-chain carnitines (B), ketone body (C), and BAs (D) in different groups. The boxes represent the middle 50% of the data, with the line inside the box indicating the median and the cross indicating the mean. The whiskers indicate the maximum and minimum values and circles indicate extreme data point. (E) Weighted correlation network analysis of metabolites. Each node represents a compound, each line represents the correlation between the two compounds, and the line weight represents the correlation coefficients

Next, we performed correlation analysis on the top 77 metabolites profiles to detect the relationship of functional compounds that changed during the process of ATB-DILI. In supplementary Fig. 8, the top 77 metabolites formed two clusters that mainly belonged to FAs metabolism and BAs metabolism and compounds of the same groups tended to aggregate with each other. We further conducted weighted co-expression network analysis to explore the potential regulatory relationship of metabolites. In Fig. 2E, the largest cluster was composed of a variety of long-chain saturated FAs, monounsaturated FAs, and polyunsaturated FAs. Within the FAs cluster, intra-correlation was remarkably high. The second cluster contained multiple long-chain carnitines and ketone body, which highly correlated with each other. Furthermore, primary and secondary BAs also positively correlated with each other. All these results suggested that metabolic pathways, including long-chain FAs synthesis, elongation, desaturation, β-oxidation, ketone body synthesis, and BAs synthesis, are up-regulated with ATB-DILI occurrence and progression.

Metabolic dynamics and pathway analysis of ATB-DILI

We longitudinally examined global dynamic changes of all 895 metabolites. There were 8 clusters observed in the ATB-DILI group and ATB-Ctrl group and metabolites of each cluster showed unique metabolic characteristics (Fig. 3A and B). Metabolites that involved in the cluster 2, 3, 5, 7, 8 in the ATB-DILI group showed significant metabolic inflection point at T4, we thus speculated that compounds of the clusters were significantly associated with the occurrence and progression of ATB-DILI. Further enrichment and pathway analysis revealed that alpha linolenic acid and linolenic acid metabolism, fatty acid synthesis, carnitines synthesis, and ketone body metabolism were markedly enriched (Fig. 3C and D). These findings suggested that lipid metabolism and energy metabolic balance might be disrupted by antituberculosis drugs. Under the sustained effect of antituberculosis drugs, BA synthesis, the down-stream metabolic pathway of lipids, could also be significantly up-regulated in the pathological state.

Fig. 3

Enrichment, and pathway analysis of the metabolites. (A-B) Cluster analysis showed the temporal changes of metabolite intensity during antituberculosis treatment in the ATB-DILI group (A) and ATB-Ctrl group (B). (C-D) Enrichment and pathway analysis of metabolites that showed significant metabolic inflection point at T4

Prediction biomarkers identification for ATB-DILI

To find candidate biomarkers of ATB-DILI at the early stage, we further applied the random forest (RF) classification algorithm to select the significant variables among the top 77 metabolites. From the longitudinal perspective, we identified metabolites that differentiated samples between multiple timepoints (T2-T4) of disease progression and baseline (T1) in the ATB-DILI group. For this analysis, the prediction accuracy of RF model was 88.89%, 83.33% and 94.44% at T2-T4, respectively. In supplementary Fig. 9A-9 C, top 30 metabolites were listed in order of their importance to the classification scheme. Among them, glycohyocholate (GHCA) and taurohyocholate (THCA) ranked within the top 4 of the lists at different timepoints (T2-T4). Therefore, we speculated that GHCA and THCA could potentially serve as biomarkers for the disease. From the horizontal perspective, we identified compounds that differentiated individuals between ATB-Ctrl group and ATB-DILI group at T2-T4. For this analysis, the prediction accuracy of RF model was 61.11%, 66.67% and 83.33% at T2-T4, respectively. In supplementary Fig. 10A-10 C, we observed that FAs play a major role at T2 and BAs play a vital role at T4. Thus, 5 FAs including myristate (14:0) (MA), palmitoleate (16:1n7) (POA), linoleate (18:2n6) (LA), arachidonate (20:4n6) (ARA), oleate/vaccenate (18:1) (OA) and 6 BAs including glycochenodeoxycholic acid (GCDCA), taurochenodeoxycholic acid (TCDCA), glycocholic acid (GCA), taurocholic acid (TCA), GHCA, and THCA were involved in the following targeted verification (supplementary Table 9).

Candidate FAs validation in longitudinal cohort and cross-sectional cohort

To validate the predictive capability and universal applicability of candidate biomarkers, we further performed quantitative detection for each biomarker in the longitudinal and cross-sectional cohorts, respectively. Demographics and clinical parameters of the independent cross-sectional cohort are displayed in the Table 2.

Table 2 Clinical features of the cross-sectional cohort

In Fig. 4A, we observed that the FAs abundance increased dramatically over time in the ATB-DILI group. To evaluate the ability of 5 FAs to recognize ATB-DILI cases at the early stage, we next performed receiver operating characteristic curve (ROC) analysis at T2-T4, respectively. The area under the curves (AUC) of 5 FAs ranged from 0.63 to 0.69 at T2 and LA had the highest AUC of 0.69; the AUC fluctuated between 0.62 and 0.77 at T3 and MA showed the best performance with AUC of 0.77; the AUC ranged from 0.65 to 0.75 at T4 and ARA had the best diagnostic efficiency with an AUC of 0.75 (Fig. 4B).

Fig. 4

Concentration of candidate biomarkers and ROC analysis in the longitudinal cohort. Plasma concentration (A) and ROC analysis (B) of 5 FAs in the longitudinal cohort. Plasma concentration (C) and ROC analysis (D) of 6 BAs in the longitudinal cohort. Data are expressed as mean with SEM (A, C). * indicates p < 0.05

We then explored the abundance changes of FAs in the independent cross-sectional cohort. In Fig. 5A, we observed that the concentrations of FAs in the CS-DILI group were significantly higher than those of the CS-Ctrl group and HC group (p < 0.05). The ROC analysis showed that the AUC of OA, POA, ARA, LA, and MA were 0.62, 0.61, 0.60, 0.59, and 0.55, respectively.

Fig. 5

Concentration of candidate biomarkers and ROC analysis in the cross-sectional cohort. Plasma concentration and ROC analysis of 5 FAs (A) and 6 BAs (B) in the cross-sectional cohort. Data are expressed as median with interquartile range, * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001, **** indicates p < 0.0001

These results further validated that the level of FAs is up-regulated during the progression of ATB-DILI. And 5 FAs had the ability to identify ATB-DILI subjects at the early stage (T2, T3) while ALT remained low level (< 3×ULN) (supplementary Fig. 11).

Candidate BAs validation in longitudinal cohort and cross-sectional cohort

In the ATB-DILI group, the abundance of GCA, TCA, GCDCA, TCDCA was higher than that in ATB-Ctrl group at T4 while the concentrations of GHCA and THCA were higher than those in ATB-Ctrl group at T2-T4 (Fig. 4C). The ROC analysis revealed that 6 BAs had the best diagnostic performance at T4 and the AUCs of GCA, TCA, GCDCA, TCDCA, GHCA, THCA were 0.73, 0.73, 0.62, 0.64, 0.72, and 0.68, respectively (Fig. 4D). The prediction capability of 6 BAs was limited at the early stage (T2, T3) and the AUCs were just in the range of 0.51∼0.59. Consistent with the observation in the longitudinal cohort, the levels of GCA, TCA, GCDCA, TCDCA, GHCA, THCA dramatically increased in CS-DILI group (Fig. 5B). The ROC analysis showed that 6 BAs had a good capability to differentiate patients with ATB-DILI from controls. The AUCs of GCA, TCA, GCDCA, TCDCA, GHCA, THCA were 0.75, 0.74, 0.70, 0.70, 0.64, and 0.65, respectively.

These results further validated that the level of 6 BAs was significantly up-regulated in ABT-DILI subjects. And 6 BAs had a better diagnostic performance when ATB-DILI occurred. We have also tried logistic regression with all or part of biomarker candidates, but the AUC is lower than a single biomarker (data not shown).