Understanding how a drug interacts and affects biological systems to induce hepatotoxicity with the drug-induced transcriptome data is a goal of this study. The perturbed biological mechanisms by dysregulated genes can be interpreted as a degree of hepatotoxicity at transcriptomic level. However, despite the abundance of drug response and gene expression data, measuring the degree of drug-induced hepatotoxicity faces three major challenges.

1. Data availability varies greatly depending on drug treatment conditions. For example, the LINCS dataset includes gene expression values for 12,328 genes across an average of 16 samples (with a standard deviation of 30) for varying dose and time point combinations (Supplementary Figure S1). It causes a high-dimensional low-sample issue and hinders accurate interpretation of the MoA of drugs.

To address these challenges, we propose MDTR, a knowledge-guided Multi-Dimensional Transcriptomic Ruler for quantifying the degree of drug-induced hepatotoxicity at the transcriptomic level. We first compiled five biological mechanisms of hepatotoxicity from the recent literature (Andrade et al., 2019; Weaver et al., 2020; Han et al., 2020). To translate transcriptome data to the five hepatotoxicity mechanisms, MDTR incorporates KEGG pathways as representative of mechanisms, mapping transcriptome data to these pathways and subsequently aggregating them for each of the five hepatotoxicity mechanisms. This involves three steps: (1) Identifying the most-perturbed transcriptomic samples as outliers using Dual-SVDD. The rationale behind this assumption is that toxic signatures are distinct from all remaining transcriptomic samples. (2) Constructing a transcriptomic embedding space using the Radial Basis Function (RBF) kernel and defining a ruler in the kernel space. (3) Extending the ruler into a five-dimensional radar chart for interpreting hepatotoxicity based on the knowledge-guided biological mechanisms.

As a result, MDTR represents the degree of DILI as the distances that measure the perturbation of the biological mechanisms under the drug treated environment. MDTR distance outperforms existing methods, which measure distances between transcriptomic samples, in reflecting the dose-dependent effects of drug on liver injury. In addition to its quantitative capabilities, MDTR provides interpretive power for understanding the MoA. We also provide a website for calculating and visualizing the five-dimension radar chart using drug-treated transcriptomic data input by scientists.

As mentioned in challenge #2, the signatures that can be utilized to measure drug-induced hepatotoxicity at the cellular level have not yet been clearly defined. To bridge the current gap in available data and knowledge, we explore hepatotoxic signatures through two approaches: (i) data-driven view: identification of potential toxic samples by Dual-SVDD. (ii) knowledge-driven view: literature mining of biological mechanisms related to drug-induced hepatotoxicity.

Through the exploration of hepatotoxic signatures from both data-driven and knowledge-driven perspectives, we identified PT samples with a high likelihood of harboring drug-induced hepatotoxic signatures, as well as biological mechanisms associated with hepatotoxicity and corresponding pathways. Leveraging this valuable information, we introduced MDTR, quantifies the degree of hepatotoxicity as dysfunctions of biological mechanisms. The MDTR is represented as a five-dimensional radar chart, with each dimension corresponding to a biological mechanism, calculated through mechanism-level toxic distances.

Formally, given a biological mechanism M = { p 1 , … , p k } of k pathways, the mechanism-level toxic distance D M ( x ) (Equation 3) is defined as a sum of multiple weighted pathway-level distances as below: D M x = ∑ p ∈ M w p , x d p x (3) where w ( p , x ) (Equation 5) is a weight of pathway p with respect to x and d p ( x ) (Equation 4) is a pathway-level distance that measures the activity of how the pathway p is dysregulated on x , which is calculated on the non-linear kernel space to learn complex interactions of genes.

Given a pathway p = { g 1 , … , g l p } consisting of l p genes and PT samples X p P T , the toxicity kernel space H p ∈ R d is derived from X p P T using KernelPCA with RBF kernel, implemented via the Python Scikit-learn package. The gamma parameter γ is set to 1 / l , while all other parameters remain at their default settings. To construct the toxicity kernel space H p , the dimensionality d of H p is determined by the first min ( 0.1 ∗ l p , 3 ) principal components. By considering the relationships between transcriptomic samples in the latent space represented by genes within pathways, we can alleviate the high dimensionality issue mentioned in challenge #1. After the construction of toxicity kernel space, PT samples X p P T and non-PT samples X p non−PT are mapped into the toxicity kernel space H p with sample distribution Q p ∈ R d . Then, a pathway-level distance d p ( x ) is calculated as the Mahalanobis distance in the toxicity kernel space of the mapped samples: d p x = z − μ p T S p − 1 z − μ p (4) where z represents the latent embedding of the sample x mapped into the kernel space H p . μ p and S p are the mean vector and the positive-definite covariance matrix derived from Q p , respectively. In other words, once the toxicity kernel space is constructed using PT samples, all samples, including both PT and non-PT samples, are mapped into this space. The Mahalanobis distance is then calculated for each sample in this toxicity kernel space. After calculating d p ( x ) for all samples, min-max normalization is performed to adjust the scale.

The weight w ( p , x ) represents the significance of dysregulation of the pathway p in the sample x against the pathway-level distances of X p non−PT . Formally, as our observation, the distribution of { d p ( x ) | x ∈ X p non−PT } follows F-distribution. Then, the weight w ( p , x ) is calculated as: w p , x = − log F d p ; d 1 , d 2 = − log I d 1 ⋅ d p / d 1 ⋅ d p + d 2 d 1 2 , d 2 2 (5) where d p denotes d p ( x ) . d 1 and d 2 are degrees of freedom that measured from { d p ( x ) | x ∈ X p non−PT } . I is the regularized incomplete beta function.

The website (http://biohealth.snu.ac.kr/software/MDTR) is designed to measure the degree of liver toxicity of drugs by analyzing their gene expression data. The website uses MDTR to show the potential liver toxicity of drugs across the five mechanisms of hepatotoxicity.

Additionally, the website leverages the LINCS dataset to access gene expression data from 514 drugs across 70 different cell lines, specifically utilizing L1000 Level 5 data, which are generated by calculating z-scores relative to the controls. This provides a comprehensive assessment of potential hepatotoxicity at the transcriptomic level. The website is implemented using the Django 3.2 framework and Bootstrap 4.6 for the front-end. The interactive radar chart, representing the indicators of liver toxicity, is generated using Chart. js. Each axis of the radar chart represents an indicator of liver toxicity. The website is tested for compatibility and functionality on various web browsers, including Chrome, Microsoft Edge, Firefox, and Safari.

Two databases, Drug Induced Liver Injury Rank (DILIrank, Chen et al. (2016)) and LiverTox (Hoofnagle, 2013), were used to gather information on drug-induced hepatotoxicity. DILIrank categorizes drugs into four classes based on their potential for causing hepatotoxicity, while LiverTox assigns likelihood scores indicating the extent of reported liver injury cases. For this study, drugs categorized as most-DILI concern in DILIrank, and A and B in LiverTox were considered hepatotoxic. Drugs categorized as no-DILI concern in DILIrank, and D and E in LiverTox were considered non-hepatotoxic. A total of 220 hepatotoxic drugs and 402 non-hepatotoxic drugs were used in the analysis.

Drug-treated samples from the LINCS database (accession number GSE92742) were obtained for investigating drug-induced hepatotoxicity at the transcriptomic level. For the gene expression matrix of these samples, we used Level 5 data, which are generated by calculating z-scores relative to the controls. Categorization of these samples into DMSO, non-toxic, and toxic groups was guided by annotations from DILIrank and LiverTox. A total of 20,529 drug-treated samples (6,405 toxic, 11,333 non-toxic and 2,791 DMSO-treated samples) were collected. From the LINCS dataset, among the 10,174 LINCS best inference genes, 4,693 genes related to drug-induced liver injury were selected based on the Comparative Toxicogenomics Database (CTD, Davis et al. (2021)) and enrichment analysis of Gene Ontology (GO) terms. Further description and discussion of drug hepatotoxicity information, transcriptome data collection, and the selection of gene sets are provided in the Supplementary Methods section.

Among the 1,896 drug-treated transcriptome samples, Dual-SVDD selected the most perturbed samples, the extreme point of MDTR. As a result, 64 potentially toxic (PT) samples were identified (Supplementary Table S6). We showed the utility of Dual-SVDD through its ability to distinguish toxic samples and the robustness of dual-boundaries (Supplementary Figure S2). In addition, to evaluate the reliability of the proposed measurement, we conducted self-validation by partitioning the total dataset and performing cross-check analyses (Supplementary Figure S3). As a result, our measurement method showed consistent distances regardless of the configuration of the data.

Based on the PT samples and the curated biological mechanisms, MDTR measures the degree of DILI. Figure 2A shows a distribution of maximum distances, representing the density distribution of maximum values on the MDTR of the samples. Notably, significant statistical differences were observed between the sample groups identified by Dual-SVDD. This indicates that the MDTR distance measured from PT samples meaningfully reflects the perturbation of gene expression affected by drug treatments. Among the drug-induced transcriptome data, Figure 2B illustrates examples of radar charts from MDTR. The red-colored radar chart (Figure 2B, right), representing the treated compound labels as structural toxicity as ‘Toxic’, shows greater distances in comparison to the other two examples. Interestingly, even when the chemical structural information is identical and labeled as ‘Non-toxic’, the MDTR results of the two drug-induced transcriptome samples (green vs. blue-colored) show obvious differences. This indicates the importance of analyzing drug hepatotoxicity at the transcriptome level.

We further investigated the capability of MDTR in identifying samples treated with hepatotoxic drugs. We screened top-rankded drugs for liver toxicity, labeled as ‘Most-DILI’ in DILIrank or assigned category ‘A’ in LiverTox, as well as drugs with an absence of reported toxicity, labeled as ‘No-DILI’ in DILIrank or category ‘E’ in LiverTox. Then, we used samples treated with the selected drugs in liver cell lines, including HEPG2, HUH7, and PHH. This resulted in the identification of 10 toxic and 43 non-toxic samples, respectively. When comparing the distances between the two sample groups, we observed larger distances for toxic samples across all mechanisms, with statistically significant differences found in all mechanisms except for Bile acids accumulation (Figure 3). Interestingly, similar results were observed when the experiment was conducted not only in liver cell lines but also in the entire cell lines belonging to the data used (Supplementary Figure S4). Therefore, we showed the effectiveness of the multi-dimensional ruler in stratifying toxic and non-toxic samples at the transcriptomic level.

Our research goal is not only to measure the five-dimensional toxic mechanisms in drug-treated transcriptomic samples, but also to facilitate the understanding of toxicity mechanisms using the multi-dimensional ruler in the form of the radar chart. We explored the drug treated environment, especially the sample distance according to the drug dose. For dose-dependence analysis, a total of 198 combinations of 33 drugs with five or more drug dose points in fixed cell line and time point were used. Under the assumption that toxicity increases with increasing drug dose (Schenker et al., 1999; Hoofnagle and Björnsson, 2019), the correlation between dose and distance was calculated, and then evaluated that higher correlation resulted in better results. We compared MDTR distance with three methods: (1) Mathematical distances (Cosine, Euclidean, and Mahalanobis) calculated from the expression values of all genes in LINCS, (2) Transcriptomic Signature Distance (TSD, Manatakis et al. (2020)), and (3) Pathway Activity Score Learning (PASL, Karagiannaki et al. (2023)) with mathematical distances similar to (1) (details in Supplementary Methods).

As shown in Table 2, MDTR outperformed the other distance methods, achieving 70% positive ratio. Additionally, MDTR consistently showed the highest positive correlation ratio even when the number of drug dose points was three or four (Supplementary Table S7). Notably, the Mahalanobis distance yielded better results than Cosine and Euclidean distances. In addition, distances calculated based on biological information (MDTR and PASL) exhibited better performances compared to using the entire genes without the prior knowledge.

The main advantages of MDTR are that it not only provides distances indicating the degree of toxicity depending on the drug-treated environment but also allows interpretation of the mechanisms of liver toxicity. Among drugs that show a positive correlation with MDTR and dose, we conducted case studies for the following three drugs: Doxorubicin, Mitoxantrone, and Rosiglitazone.

Case 1

We investigated Doxorubicin, a chemotherapy medication widely used for the cancer treatments, including breast and bladder cancer, while also known for its hepatotoxicity (Prasanna et al., 2020). In the MCF7 cell line and at the 6-h time point, Doxorubicin exhibited dose-dependent toxicity signature (Figure 4, ρ = 0.93, p = 1.1e-6). These positive correlations were observed in 9 out of the 10 Doxorubicin-treated combinations (Supplementary Table S8). In particular, the radar chart revealed distinct changes across different dose points in two mechanisms: Mitochondrial dysfunction and Bile acids accumulation. These findings support the understanding of the potential association between bile acid metabolism and Doxorubicin sensitivity (Chen et al., 2017), as well as the preferential accumulation of Doxorubicin in the mitochondria and nucleus (Wallace et al., 2020).

Case 2

We investigated the potential toxicity of Mitoxantrone, a cytotoxic agent and anthracenedione analogue of Doxorubicin (Fox, 2004), which exhibits similar metabolic pathways (Paciucci and Sklarin, 1986; Llesuy and Arnaiz, 1990; Rossato et al., 2014a). Mitoxantrone, despite being classified with a likelihood score of D in the LiverTox database, showed a positive correlation between distance and dose across all combinations of MCF7, PC3, and A549 cell lines and 6 and 24 h time points, similar to Doxorubicin (Supplementary Table S9). For example, Figure 4B illustrates the dose-dependent changes in distance for Mitoxantrone in the PC3 cell line at the 24-hour time point ( ρ = 0.87, p = 0.012).

In particular, the toxicity of the sample showed a non-monotonic tendency that decreased in the early stage of dose escalation (from 0.1 to 0.5 μ M) and then increased again after the 1uM dose (Tang et al., 2023). These results showed that MDTR is capable to capture the pattern of drugs showing such hormesis. Overall dose-distance relationship suggest that Mitoxantrone not only induces functional impairments observed with Doxorubicin, such as bile acids accumulation and mitochondrial dysfunction, but also affects lipid metabolism (Kharasch and Novak, 1985; Rossato et al., 2014b).

Furthermore, we performed the MDTR on an unseen drug, Rosiglitazone, a drug that has a lack of medical reports concerning liver toxicity and was not utilized in the generation of dual-boundaries or the calculation of the ruler.

Case 3

We investigated the liver toxicity of Rosiglitazone, a thiazolidinedione drug used for the treatment of type 2 diabetes treatment but discontinued in many countries due to cardiovascular risks (Nissen and Wolski, 2007). While its counterpart Troglitazone has been withdrawn from the market due to severe hepatotoxicity (Lloyd et al., 2002; Kaplowitz, 2005), limited research exists on the hepatotoxic effects of Rosiglitazone. To address this gap, we employed the multi-dimensional ruler to evaluate its toxicity. In the PC3 cell line at the 24-h time point, Figure 4C demonstrates an increase in toxic mechanisms corresponding to higher drug doses. Specifically, significant variations were observed in Bile acids accumulation and Lipid metabolism in response to varying doses. These findings align with known issues associated with Rosiglitazone, such as disturbances in fatty acid and triglyceride metabolism (Gershell, 2005; Tan et al., 2005; Watkins et al., 2002), as well as intrahepatic cholestasis (Zhang et al., 2020).

Through comparative and case study results, we show that MDTR is an interpretable method that can capture both known and potentially mechanisms of toxicity, while exhibiting a dose-dependent distance to hepatotoxicity.

Oxidative stress, one of the axis of MDTR, arises from an imbalance between oxygen-reactive species (ROS) generation and accumulation (Figure 5A, Jaeschke et al. (2002), Jaeschke et al. (2012)). The accumulation of ROS, triggered by factors such as AMP-activated protein kinase (AMPK) activation (Kang et al., 2016; Steinberg and Hardie, 2022) and an imbalance in glutathione (GSH) and glutathione disulfide (GSSG) levels (DeLeve and Kaplowitz, 1991; Yuan and Kaplowitz, 2009), promotes inflammation, disrupts cellular energy regulation, depletes antioxidants, impairs protein handling, and induces hepatocellular apoptosis, ultimately leading to liver damage.

To explore the biological implications of pathway-level distances within the mechanism, we focused on the Chemical carcinogenesis - ROS pathway (hsa05208, hereafter referred to as the ROS pathway), which is one of the six pathways chosen for the oxidative stress mechanism. We divided the samples evenly into four groups (Low, Middle1, Middle2, and High) based on each pathway distance. Then, we performed co-expression network analyses based on the expression values of samples within each group (Figure 5B). As a result, larger distances correspond to large co-expression network sizes, indicating that increased distance is associated with more interactions between genes and, simultaneously, a higher degree of gene dysregulation. For example, during the transition from the Low group to the Middle1 group, downregulation of MAPK8 (also known as JNK) may hinder the expression of the AP-1 transcription factor, potentially impacting the MAPK signaling pathway (Turner et al., 2014). Furthermore, as the network progressed from the Middle1 to Middle2 and High groups, upregulation of genes associated with PI3K signaling (e.g., PIK3CA or PIK3CD) and NADPH oxidase (e.g., NOX4) was observed, indicating ROS production (Koundouros and Poulogiannis, 2018). These findings validate the efficacy of the proposed pathway-level distance in capturing biological activity. Similar results were observed for the remaining pathways (Supplementary Figures S5-S9).

We next calculated the mechanism-level distance of oxidative stress as a weighted sum of the six related pathways, including the ROS pathway. To explore the association between oxidative stress signaling and the mechanism-level distance of samples, we performed GO enrichment analysis using commonly perturbed genes from the respective group. Figure 5C and Supplementary Figure S10 illustrate that as the distance increases, the GO terms related to oxidative stress, such as response to endoplasmic reticulum stress and cellular response to oxidative stress, exhibit more significant and abundant enrichment. These GO enrichment analysis results were consistently observed across all the mechanisms (Supplementary Figures S11-S14). Thus, our findings demonstrate that both pathway- and mechanism-level distance effectively reflect the activity level of the corresponding mechanism as the distance increases.

Our model is implemented as a freely available website (http://biohealth.snu.ac.kr/software/MDTR). The website allows for our MDTR to be utilized from two perspectives (Supplementary Figure S15): (1) detection of hepatotoxicity from the user’s drug-induced transcriptome data, (2) investigation of hepatotoxicity through drug-induced transcriptome data in the LINCS database by selecting any combinations of drugs, tissues, and cell lines.