This prospective cohort study included patients with tuberculosis who visited the tuberculosis clinic of West China Hospital of Sichuan University from March 2021 to December 2021. The study was approved by the Ethics Committee of West China Hospital of Sichuan University. All research subjects were required to sign a written informed consent form by themselves or their representatives before being included in the study. Demographic datasets of patients with laboratory test data were obtained through electronic medical records and questionnaires.

Inclusion criteria are as follows: 1) age≥16 years and <80 years old; 2) newly diagnosed TB patients, including etiologically confirmed, pathologically confirmed and clinically diagnosed cases; 3) standard first-line anti-TB treatment regimens (including 2-month HRZE intensive treatment and at least 4 months of HRE consolidation therapy), and can be followed up regularly; 4) Han nationality in Southwest China; 5) voluntarily participate in this study and sign the informed consent form when included in the study. Those who do not meet the above diagnostic criteria are referred to as non-ATB-DILI.

The exclusion criteria were as follows: 1) abnormal liver function at baseline; 2) concomitant liver diseases such as alcoholic hepatitis, viral hepatitis or liver cirrhosis; 3) taking immunosuppressive drugs, antitumor drugs, and acetaminophen and other drugs that may cause liver damage; 4) patients with diabetes, autoimmune diseases, malignant tumors, or tuberculosis with severe heart, lung, and renal insufficiency; and 5) patients with HIV infection or who died during follow-up from causes other than adverse drug reactions.

The diagnostic criteria for ATB-DILI used in this study are (9, 34–36) as follows: alanine aminotransferase (ALT)≥3 normal upper limit of normal value (ULN) and/or total bilirubin (TBil)≥ 2ULN; or aspartate aminotransferase (AST) or alkaline phosphatase (ALP) and TBil are elevated at the same time, and at least one of them is ≥2 ULN. In addition, considering that the updated Roussel Uclaf Causality Assessment Method (RUCAM) was the recognized standard for the diagnosis of DILI, no matter what drug was used (37). Here, for patients who met the DILI diagnostic criteria (ALT≥5ULN or ALP≥2ULN) in the updated RUCAM, we conducted a subgroup analysis. All patients diagnosed with ATB-DILI completed a causality assessment using the updated RUCAM scale (37). Only patients with a RUCAM score of 6 or more were accepted in this study.

Urine nontargeted metabolomics and microbiome analyses were performed in ATB-DILI and gender and age-matched non-DILI patients, and patients with urinary system diseases were excluded. Urine samples were self-collected from the subjects according to the provided instructions prior to antituberculosis treatment and immediately sent to the laboratory. Once in the laboratory, samples were stored at −80°C until metabolomic and microbiome analysis.

Urine is the common sample type used to perform metabolomics studies (38, 39). Compared to other samples, urine has easy sampling, low protein levels and less complexity (38). Also the urine metabolites are products of normal and abnormal cellular biological processes and can reflect a wide range of phenotypes including genetic modifications (38). Therefore, urine is more advantageous compared to other sample types and was used as a study sample in this study.

Clean midstream urine from patients before medication was collected, divided into three 1 ml aliquots, and immediately stored at −80°C. We discarded samples that were at room temperature for >2 hours.

Metabolite extraction was primarily performed according to previously reported methods (40, 41). In short, 100 µL samples were extracted by directly adding 300 µL of precooled methanol and acetonitrile (2:1, v/v), and internal standards mix (contains: L-Leucine-d3, L-Phenylalanine (13C9, 99%), L-Tryptophan-d5, Progesterone-2,3,4-13C3) were added for quality control of sample preparation. After vortexing for 1 min and incubating at -20°C for 2 h, the samples were centrifuged for 20 min at 4000 rpm, and the supernatant was then transferred for vacuum freeze drying. The metabolites were resuspended in 150 µL of 50% methanol and centrifuged for 30 min at 4000 rpm, and the supernatants were transferred to autosampler vials for LC‒MS analysis. A quality control (QC) sample was prepared by pooling the same volume of each sample to evaluate the reproducibility of the whole LC‒MS analysis.

This experiment used a Waters 2D UPLC (Waters, USA) tandem Q Exactive HF high resolution mass spectrometer (Thermo Fisher Scientific, USA) for separation and detection of metabolites. To provide more reliable experimental results during instrument testing, the samples are randomly ordered to reduce system errors. A QC sample was interspersed for every 10 samples.

The samples were analyzed on a Waters 2D UPLC (Waters, USA), coupled to a Q-Exactive mass spectrometer (Thermo Fisher Scientific, USA) with a heated electrospray ionization (HESI) source and controlled by the Xcalibur 2.3 software program (Thermo Fisher Scientific, Waltham, MA, USA). Chromatographic separation was performed on a Waters ACQUITY UPLC BEH C18 column (1.7 μm, 2.1 mm × 100 mm, Waters, USA), and the column temperature was maintained at 45°C. The mobile phase consisted of 0.1% formic acid (A) and acetonitrile (B) in the positive mode, and in the negative mode, the mobile phase consisted of 10 mM ammonium formate (A) and acetonitrile (B). The gradient conditions were as follows: 0-1 min, 2% B; 1-9 min, 2%-98% B; 9-12 min, 98% B; 12-12.1 min, 98% B to 2% B; and 12.1-15 min, 2% B. The flow rate was 0.35 mL/min and the injection volume was 5 μL.

The mass spectrometric settings for positive/negative ionization modes (ESI+/-) were as follows: spray voltage, 3.8/−3.2 kV; sheath gas flow rate, 40 arbitrary units (arb); aux gas flow rate, 10 arb; aux gas heater temperature, 350°C; capillary temperature, 320°C. The full scan range was 70–1050 m/z with a resolution of 70000, and the automatic gain control (AGC) target for MS acquisitions was set to 3e6 with a maximum ion injection time of 100 ms. The top 3 precursors were selected for subsequent MSMS fragmentation with a maximum ion injection time of 50 ms and resolution of 30,000, and the AGC was 1e5. The stepped normalized collision energy was set to 20, 40 and 60 eV.

The original data (raw file) collected byLC‒MS/MS were imported into Compound Discoverer 3.1 (Thermo Fisher Scientific, USA) for data processing, including peak extraction, retention time correction, background peak labeling, and metabolite identification. We calculate the coefficient of variation of the relative peak area in all QC samples, and delete the compounds with coefficient of variation greater than 30%. The identification of metabolites was a combined result of the BGI Metabolome Database (BMDB), mzCloud and ChemSpider (Human Metabolome Database (HMDB), Kyoto Encyclopedia of Genes and Genomes (KEGG), LipidMaps) databases. Main parameters of metabolite identification: Precursor Mass Tolerance <5 ppm, Fragment Mass Tolerance <10 ppm, RT Tolerance <0.2 min. The identification level of metabolites was divided into five confidence levels, and the credibility of Level 1 to Level 5 decreased in order. The original data exported by LipidSearch were imported into metaX for data preprocessing and subsequent analysis (42). Multivariate statistical analysis (principal component analysis (PCA) and partial least squares-discriminant analysis (PLS-DA)), and univariate analysis (fold-change, FC and Student’s t test) were combined to screen for differential metabolites between groups. Differential metabolite screening conditions: 1) variable projected importance (VIP) ≥ 1, 2) fold-change ≥ 1.2 or ≤ 0.83, 3) p-value <0.05. Metabolic pathway enrichment analysis of differential metabolites was performed based on the KEGG database.

Microbial genomic DNA extraction was performed as described previously (43). Urine microbial DNA was extracted using a Qiagen Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. Primers targeting the hypervariable V3+V4 region of the 16S gene were used to amplify the extracted DNA samples (the forward primer was 5ʹ- ACTCCTACGGGAGGCAGCA -3ʹ, and the reverse primer was 5ʹ- GGACTACHVGGGTWTCTAAT -3ʹ). All samples were sequenced via Illumina HiSeq 2500.

Cutadapt v2.6 software was used to process the raw data to obtain fragments of the target region. FLASH (Fast Length Adjustment of Short reads, v1.2.11) was used for sequence splicing, and UCHIME (v4.2.40) software was used to remove chimeras. Sequences were clustered with a 97% similarity level by using USEARCH (v7.0.1090_i86linux32) to cluster the spliced tags into OTUs. The OTU representative sequences were aligned with the database for species annotation by RDP classifier (1.9.1) software, and the confidence threshold was set to 0.6. The VennDiagram package of R (v3.1.1) software was used to display the number of common and unique OTUs for each group. Principal coordinate analysis was performed using QIIME (v1.80) software to present similarities or differences in data. Line Discriminant Analysis (LDA) Effect Size (LEFSE) was used to calculate the differences in species abundance between the two groups and then to research the biomarkers related to ATB-DILI.

Differences between two groups were compared by using Student’s t test for normal continuous variables and χ2-test for categorical variables. Differences with a p value <0.05 (two-sided) were considered statistically significant. Statistical analyses were performed using SPSS V.21.0 for Windows (SPSS, Chicago, Illinois, USA). Moreover, correlations between the microbiota and metabolites and between the metabolites and clinical parameters were analyzed. Integrating electronic medical records, metabolomic and microbiome data, and machine learning methods was used to establish a clinical prediction model of ATLI. We used four machine learning algorithms: random forest, artificial neural network, support vector machine (SVM) with the linear kernel (SVM-linear), and SVM with radial basis function kernel (SVM-rbf) (44). The stratified sampling method was used to divide the training set (80%) and the validation set (20%), and the R4.1.2 software (R Foundation for Statistical Computing, Vienna, Austria) was used for data screening and model building. The importance of each feature in the occurrence of ATB-DILI was scored, and area under receiver operating characteristic (ROC) curves were employed to assess the accuracy of the models.

A total of 74 patients diagnosed with TB were recruited for this study from March 2021 to December 2021 at West China Hospital of Sichuan University (Sichuan Province, China). Finally, 5 patients were lost to follow-up. Of the remaining 69 patients, 16 (23.19%) developed ATB-DILI after antituberculosis treatment ( Supplementary Figure 1 ). The general clinical characteristics of the two groups and the results of liver function tests when DILI occurred are shown in Table 1 . Compared with the non-ATB-DILI group, the levels of albumin (43.2(40.1-44.4)g/L vs. 44.9(42.2-46.8)g/L, P: 0.033) and hemoglobin (125.5(118.5-135.8)g/L vs. 137.0(128.5-145.5)g/L, P: 0.019) were significantly lower in the ATB-DILI group. No significant differences were observed in other baseline characteristics between the two groups of patients (P>0.05). The median time to DILI occurred on day 29 after taking anti-TB drugs ( Table 1 ).

After excluding patients with urinary system diseases, 14 ATB-DILI patients and 30 age- and sex- matched non-ATB-DILI patients were included for urinary metabolomics and microbiome analysis ( Supplementary Figure 1 ). It is important to note that of these 14 patients, 8 met the definition of liver adaptation (45). The other 6 patients with ALT≥5 times the upper limit of normal were stopped using antituberculosis drugs according to Chinese guidelines, so it was hard to distinguish which were liver adaptation (36).

As shown in Table 2 , there were no significant differences in sex, age, body weight, BMI, body mass index (BMI), smoking, drinking or tuberculosis site between the two groups of patients who underwent urine nontargeted metabolome and microbiome analysis (P>0.05). All participants had normal liver function before anti-tuberculosis drug ingestion. It was suggested that the general conditions of the two groups were consistent and comparable.

PCA and OPLS-DA were performed for both positive ion mode (ESI+) and negative ion mode (ESI−). As shown in the figure ( Figures 1A, B ), the QC samples (blue circles) were significantly aggregated, indicating that the instrument was stable and that the reproducibility of the acquired data was good. The ATB-DILI (n=14, red circles) and tolerance groups (n=30, green circles) were not well separated in PCA. As shown ( Figures 1C, D ), the PLS-DA model clearly separated the ATB-DILI and non-ATB-DILI groups in both ionization modes. Differential metabolites between the two groups were screened according to multivariate and univariate statistical significance criteria (VIP≥1, FC≥1.2 or ≤ 0.83, and P<0.05). In general, there were 1256 urine differential metabolites screened in the positive ion mode and 334 in the negative ion mode ( Figure 2 ). Finally, 28 differential metabolites with secondary classification names and reliable identification results (Level 1-3) were selected ( Table 3 ), including choline, cherry base, N-acetyl, pseudohadine, N8-acetyl spermamine, glycolic acid, etc. As shown in Table 4 , a total of 7 significant enrichment pathways for differential metabolites were found in both positive and negative ion modes. The differential metabolites were mainly involved in the metabolism of bile secretion, nicotinate and nicotinamide metabolism, tryptophan metabolism, ABC transporters, neuroactive ligand‒receptor interaction, arginine and proline metabolism, and porphyrin and chlorophyll metabolism (P<0.05, Count≥2) ( Table 4 ).

Correlation analysis was conducted between urine differential metabolites and clinical data, including baseline ALT, AST, TBIL, Alkaline phosphatase, hemoglobin, uric acid and albumin. We found that many different metabolites were significantly correlated with clinical data ( Supplementary Table 1 ). The urine differential metabolite 11 dehydrothromboxane B2 was positively correlated with the baseline total bilirubin concentration, while the urine differential metabolite N8-acetylspermidine was negatively correlated with the hemoglobin content, and uric acid was also negatively correlated with the baseline serum uric acid level ( Supplementary Table 2 ).

As shown in the Figure ( Figure 3A ), 1079 OTUs were shared between the ATB-DILI group and the non-ATB-DILI group, 607 OTUs were unique to the non-ATB-DILI group, and the other 189 OTUs were unique to the ATB-DILI group. The Shannon curve ( Figure 3B ) shows that the amount of sequencing data in this study was large enough to reflect the vast majority of microbial information in the sample. The top 10 key species between the two groups are shown in Figure 3C . Weighted UniFrac principal coordinate analysis (PCoA) was applied to detect the changes in microbial community structures ( Supplementary Figure 3 ). The results indicate that the ATB-DILI group and the control group were significantly separated along the PC2 axis, which explained 19.79% of the total variation. LEFSE analysis was used to determine the key attribute differences between the two groups. The differential microbiota (LDA score>2) screened between the two groups were Negativicoccus and Actinotignum, which were all upregulated in the ATB-DILI group ( Figure 3D ).

We further investigated the correlation of urinary differential metabolites with altered urinary microbiota. Significant correlations were found between some differential metabolites and microbial groups by calculating rank correlation coefficients ( Figure 4 and Supplementary Table 2 ). Carbendazim was positively correlated with synergistia but negatively correlated with mollicutes (p<0.05) ( Supplementary Table 2 ). D-(-)-lyxose was positively correlated with four microbial groups, including synergistia, ktedonobacteria, fibrobacteria and fusobacteriia (p<0.05) ( Supplementary Table 2 ). Altogether, these results showed that distinctive metabolites were closely related to urinary microbiome variation, and distinctive metabolites and microbiomes were closely related to the occurrence of ATB-DILI.

According to RUCAM criteria, metabolome and microbiome analysis were performed between the ALT≥5 ULN group and normal patients. Significant metabolic differences were observed between the two groups, there were 1122 different metabolites were screened in positive ion mode and 386 in negative ion mode ( Supplementary Figure 3 ). Finally, 26 different metabolites were selected, including choline, 11-dehydrothromboxane b2, and N8-acetylspermidine. ( Supplementary Table 3 ). Consistent with our results in Section 3.2, 8 common different metabolites were found in the subgroup analysis to be related to liver injury after medication, especially when ALT>5ULN occurred ( Table 3 and Supplementary Table 3 ). The eight differential metabolites were choline, N8-acetylspermidine, carbendazim, N-acetylputrescine, 1-methylnicotinamide, creatine, porphobilinogen, and nonanoic acid ( Table 3 and Supplementary Table 3 ). And these different metabolites had the same label direction in the two groups of patients with liver injury ( Table 3 and Supplementary Table 3 ).

There were 795 OTUs shared between the DILI patients with ALT ≥ 5 ULN and normal patients, 68 OTUs were unique to the DILI group, and the other 891 OTUs were unique to the control group ( Supplementary Figure 4-A ). The top 10 key species between the two groups are shown in Supplementary Figure 4-B . Finally, 3 differential microbiotas (LDA score>2) were found between the two groups ( Supplementary Figure 4-B ). The Actinotignum was down regulated in DILI group, while the Bradyrhizobiaceae, and Bradyrhizobium were upregulated in the non-DILI group ( Supplementary Figures 4-C ). Combined with the results in Section 3.4, we have sufficient evidence to show that Actinotignum was closely related to the occurrence of liver injury after medication, regardless of the DILI standards.

Random forest analysis was performed on the screened differential metabolites ( Table 3 ), differential microbiota ( Figure 2 ), and relevant clinical data of 44 patients. For clinical characteristics, we included albumin and hemoglobin, which were significantly different between the two groups, as well as other factors that may be associated with the occurrence of ATB-DILI (including age, sex, BMI, baseline ALT, AST, and TBil). A total of 38 variables are included. When ntree=500 and mtry=6, the model reaches the optimum. The score of the 38 variables was shown in Figure 5A . The larger the absolute value is, the greater the importance of the indicator. After sorting the variables from high to low according to the absolute value, the cross-validation curve was obtained by performing tenfold cross-validation repeated 5 times ( Supplementary Figure 5 ). The top 10 variables were selected for model building with the lowest error ( Supplementary Figure 5 ). The area under the ROC curve of the four models were shown in Table 5 . At training set, the random forest model performed significantly better than the remaining three models (area under the curve 0.98 vs. 0.87 (ANN), 0.89 (SVM-linear) and 0.89 (SVM-rbf) ( Table 5 ). Overall, random forest model, artificial neural network model and two support vector machine models (both SVM-linear and SVM-rbf) all have excellent prediction value for the validation set ( Figure 5B and Table 5 ). The consistent results between the training set and the validation set indicate that those models have high accuracy for predicting the occurrence of ATB-DILI.