A549 (Cat. No. CCL-185) and HepG2 (Cat. No. HB-8065) cells were obtained from the American Type Culture Collection (Manassas, VA, USA); HT-29 (Cat. No. 1101HUM-PUMC000109) and HEK293 (Cat. No. 1101HUM-PUMC000010) cells were obtained from the Cell Resource Center of Peking Union Medical College (Beijing, China). A549 and HepG2 cells were cultured in Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 IU/mL penicillin and 100 μg/mL streptomycin; HEK293 cells were cultured in Minimum Essential Medium supplemented with 10% FBS, 100 IU/mL penicillin and 100 μg/mL streptomycin; HT-29 cells were cultured in DMEM/F-12 supplemented with 5% FBS, 100 IU/mL penicillin and 100 μg/mL streptomycin. DMEM (Cat. No. 11965500BT), Minimum Essential Medium (Cat. No. C11095500BT), FBS (Cat. No. 10099141C), penicillin streptomycin (Cat. No. 15140122), and DMEM/F-12 (Cat. No. 11320033) were purchased from Gibco (Waltham, MA, USA).

A total of 1243 compounds (purity >95%) were obtained from the compound libraries (Cat. Nos. L1000, L3400, L6000, and L4000), and the powders of the compounds, cyclosporine A (Cat. No. T0945; purity = 98%) and triptolide (Cat. No. T2179; purity = 98%), were purchased from TargetMol (Boston, MA, USA). The structure and purity of all the compounds were confirmed by nuclear magnetic resonance and high-performance liquid chromatography. The compounds used in this study were dissolved in dimethyl sulfoxide (DMSO, Cat. No. 34943; Sigma–Aldrich, St. Louis, MO, USA) and stored at −20 ℃.

A549, HEK293, HepG2, and HT-29 cells were seeded in 96-well plates and treated with the test compounds for 48 h. Cell viability was determined via the CellTiter-Glo Assay (Cat. No. G7571; Promega, Madison, WI, USA) by measuring the relative light units (RLUs) according to the manufacturer's protocol. Cells treated with DMSO (1‰, v/v) served as the vehicle control. Cell viability was calculated using Eq. (1)¹⁹:

All the experiments were performed in duplicate, and the mean values subsequently were utilized for classifier label determination.

A549, HepG2, and HT-29 cells were seeded in 6-well plates and incubated with the test compounds or DMSO (1‰, v/v). Briefly, HT-29 cells were treated with cyclosporine A at a final concentration of 10 μmol/L; A549 and HepG2 cells were treated with triptolide at final concentrations of 30 and 10 nmol/L, respectively. Total RNA was extracted using TRIzol reagent (Cat. No. 15596026; Gibco, Waltham, MA, USA) according to the manufacturer's protocol. RNA quality assessment, cDNA library construction and RNA sequencing were conducted by OE Biotech Co., Ltd. (Shanghai, China). Briefly, RNA integrity was assessed by an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The cDNA library was constructed using a VAHTS Universal V6 RNA-seq Library Prep Kit (Vazyme, Nanjing, China) according to the manufacturer's instructions and sequenced on Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA). The raw reads were first processed using fastp²⁰ to obtain the clean reads, which were subsequently mapped to the reference genome using HISAT2²¹. The FPKM²² value of each gene was calculated, and the read counts of each gene were obtained via HTSeq-count²³. Differential expression analysis was performed using DESeq2²⁴, and the significant differential expressed genes (DEGs) were identified by using the following thresholds: q value < 0.05 and fold change >2 or < 0.5.

The transcriptome data utilized in this study were downloaded from the Library of Integrated Network-based Cellular Signatures (LINCS, https://clue.io/data/CMap2020#LINCS2020)²⁵ and “sig_id” information of the files is provided in Supporting Information Table S1. The level 3 and level 5 matrixes of 1243 compounds were extracted using the parse function of the cmapPy python package²⁶. Subsequently, the signatures obtained under the same conditions were averaged using the Moderated z score (MODZ) procedure²⁵ as follows: the pairwise Spearman correlation values of the profiles were calculated for the weight of each signature, and the values with negative correlations were set as 0.01; then the weighted average of each signature was subsequently calculated for multiple expression profiles to obtain the profile for modeling. To obtain a high-level view of the distribution of the compounds in modeling data, t-SNE analysis²⁷ was performed on the 1243 signatures for non-linear dimensionality reduction and visualization.

According to the A549 cell viability, the 1243 compounds were divided into three groups according to the following ranges: less than 50%, 50%–80%, and greater than 80%. The training set and test set of this study were generated by randomly stratifying the data at a ratio of 8:2 by using the train_test_split function.

Four machine learning algorithms, namely, support vector machine (SVM), random forest (RF), extreme gradient boosting (XGBoost), and light gradient boosting machine (LightGBM), were used to train the cell viability prediction models, including the classification and the regression models. The classification models were developed by using the svm.SVC function, the RandomForestClassifier function, the xgboost python package and the lightgbm python package from the scikit-learn library, while regression models were developed by using the svm.SVR function, the RandomForestRegressor function, the xgboost python package, and the lightgbm python package. In addition, voting and stacking were applied to fuse the four classification models or regression models, as the voting model was developed by a soft voting strategy (votingClassifier), and the stacking model (StackingClassifier) was developed by using the four models as base learners, and the logistic regression model and LinearRegression model were the meta-learners for the classification models and regression models, respectively. The hyperparameter tuning of each model was developed by using the bayes_opt python package in 5-fold cross-validation. Additionally, deep learning algorithms were also used to train cell viability prediction models: the multilayer perceptron and convolutional neural network were developed by using the Keras package with the backend of TensorFlow³². The TabResnet model and TabNet model were developed by using the pytorch_widedeep package³³ and pytorch_tabnet package³⁴, respectively.

The metrics, including accuracy, precision, recall, specificity, and area under the receiver operating characteristic curve (AUROC), were calculated by sklearn.metrics and used to evaluate the performance of the classifiers. The metrics were defined using Eqs. (3)–(6):where TP, TN, FP, and FN represent the number of true positives, true negatives, false positives, and false negatives, respectively.

SHapley Additive exPlanations (SHAP) analysis was used for interpretation analysis³⁵. The SHAP value of each gene in individual samples in the training set was calculated by the KernelShap function, and the global feature explanation for the model were achieved as the mean absolute SHAP value of each gene. For interpretation of individual sample predictions, local explanations were conducted by using the shap.summary_plot function for quantification of the local contribution of each gene.

The feature importance score of an individual feature on the 80% cell viability prediction model was represented by its mean absolute SHAP value across all samples in the training set. The features with feature importance score greater than the average feature importance score were acquired. Subsequently, Pearson correlation (ρ) was calculated by using the pandas.corr function between each important feature gene and other variable genes to obtain a comprehensive gene set associated with cell viability (ρ > 0.5), which was defined as the CTS genes in this study.

The intersection of the DEGs (or the DEGs with CTS genes removal) of the compound with the genes detected in the query platform was obtained. The top 150 up- and downregulated genes were selected according to the log₂ (fold change) and used for MOA prediction via the Query platform (https://clue.io/query)²⁵. The top 10 MOA terms of the corresponding cell line based on the norm connectivity score were obtained with the gene sets before and after the removal of the CTS genes used as inputs, respectively.

Pathway enrichment analysis was performed by using the Metascape platform (http://metascape.org)³⁶ based on the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes pathways, and reactome databases, with the following thresholds: P < 0.01, a minimum count of 3, and an enrichment factor >1.5.

All data needed to evaluate the conclusions of the paper are present in the paper and/or Supporting Information The RNA sequencing (RNA-seq) dataset generated in this manuscript has been deposited in the NCBI Gene Expression Omnibus database (GSE252529). All code for model training associated with the current submission is available at https://github.com/WuYou-imm/Interpretable-cytotoxicity-prediction-model.

Cell viability prediction models were constructed in this study (Fig. 1) by analyzing a total of 1243 bioactive samples with diverse MOAs that act on 431 targets and all samples were well-distributed across the training set and test set (Supporting Information Table S1 and Fig. S1A). The transcriptomic data of A549 cells treated with the compounds and the cell viabilities under the same conditions were correlated. The transcriptomic data were obtained by screening level 5 of the LINCS database, and the multiple expression profiles of the corresponding compounds were combined by using the MODZ (Fig. S1B). The cell viability data (Table S1 and Fig. S1C) of A549 cells treated with the test compounds were determined under the corresponding conditions.

Two thresholds, 50% and 80% cell viability, were applied in this study to construct classifiers. The 50% cell viability was commonly used as the cytotoxicity threshold, and the corresponding concentration, the 50% cytotoxic concentration, was often used to represent or compare the toxicity of compounds. In addition, given that compounds with weak toxicity have limitations in drug development, 80% cell viability was also used in this study as the threshold for distinguishing weakly toxic substances. Among the 1243 compounds, 991 compounds with cell viability greater than 50% (Fig. S1D) and 781 compounds with cell viability greater than 80% (Fig. S1E) were defined as positive samples for 50% and 80% cell viability classifiers, respectively. The 1243 compounds were subsequently divided into three groups according to cell viability, with percentages ranging from less than 50%, 50%–80%, and greater than 80%. The training set (n = 994) and test set (n = 249) of this study were acquired by random stratification at a ratio of 8:2 (Fig. S1C).

Given the high dimensionality and redundancy of the transcriptome data with 12,328 gene expression values, an ensemble feature selection pipeline was used to obtain the Ensemble Feature Set (Fig. 2A). First, the MFC and SD of all genes among the samples in the training set were calculated to eliminate genes with invariable expression levels, and a total of 5135 variable genes (SD > 1 or MFC >2; Supporting Information Table S2) were obtained for subsequent analysis. Then, WGCNA, mRMR, and GSEA were applied to 5135 variable genes. By using WGCNA (Supporting Information Fig. S2), the expressed gene modules and the correlation coefficients between the modules and cell viability were quantified. The module with the highest positive coefficient of 0.58 (P = 1 × 10⁻⁹¹; the blue module in Fig. S2C) and the module with the highest negative coefficient of −0.36 (P = 3 × 10⁻³²; the turquoise module in Fig. S2C) were obtained, and the genes in these two modules formed Gene set 1 (n = 3430) (Supporting Information Table S3). The mRMR algorithm was used to obtain the feature set with the greatest correlation with the label and minimum redundancy across features by calculating the mutual information score between features and the label, as well as the redundancy of features. The features of 50% (n = 1432) (Supporting Information Table S4) and 80% cell viability (n = 636) (Supporting Information Table S5) with a score >0 serving as the threshold were obtained for Gene set 2. In addition, the GSEA method was applied based on the BP terms in the GO database to obtain Gene set 3 by analyzing the frequency of genes associated with toxic compounds (i.e., cell viability less than 50% or 80%) since more frequently, genes enriched in compounds associated with toxicity are more likely to be highly correlated with cell viability. In Gene set 3, 50% (n = 4736) (Supporting Information Table S6) and 80% cell viability (n = 4747) (Supporting Information Table S7) were obtained with C > 0 as the threshold. The RRA method was subsequently used to aggregate the three Gene sets to obtain Ensemble Feature Sets of 50% and 80% cell viability, which included 292 (Supporting Information Table S8) and 226 genes (Supporting Information Table S9), respectively. Among the features, 197 genes overlapped (Fig. 2B), implying that the two Ensemble Feature Sets were mostly involved in the cell cycle (dots in Steel Blue, Fig. 2C and D) and protein metabolism (dots in Light Sea Green, Fig. 2C and D); moreover, the difference between the two Ensemble Feature Sets was mainly reflected in the presence of cell death pathways associated with 50% cell viability, and cell proliferation genes were associated with 80% cell viability (Fig. 2C and D). Toxic compounds could suppress cell proliferation, affect metabolism, and/or induce cell death. In the event of weak toxicity, the cellular oxidative stress response is triggered as a means to counteract cellular oxidative imbalance, which is largely reversible and serves as an adaptive stress response^37—39. However, severe toxic exposure significantly affects metabolic function and leads to heightened cell damage, inflammation, and cell death^37—39. According to the results of the enrichment analysis of the two Ensemble Feature Sets, the majority of the pathways were in common, but some specific pathways were different due to the distinct toxicity events of the compounds, which had varying levels of toxicity.

To construct models with a cell viability threshold of 50% or 80%, a sample with a cell viability greater than the threshold was considered a positive sample. Four models, SVM, RF, XGBoost, and LightGBM, were trained via hyperparameter optimization through Bayesian optimization via 5-fold cross-validation (Supporting Information Tables S10 and S11); subsequently, two ensemble models, voting and stacking, were constructed by combining the four models. Finally, the performance of the six prediction models for each threshold was evaluated on the test set, with the AUROC score serving as the primary metric. The results showed that both the 50% and 80% cell viability prediction models exhibited superior performance in predicting positive samples, with AUROC values ranged from 0.8823 to 0.9011 (Supporting Information Table S12) and 0.8198 to 0.8408 (Supporting Information Table S13), respectively (Fig. 3). The 50% cell viability is the most common indicator for comparing differences in compound toxicity or tumor sensitivity. Among the 50% cell viability prediction models, the voting model had the highest AUROC value (0.9011) and was regarded as the optimal model for 50% cell viability prediction (Fig. 3A). The accuracy, precision, recall, and specificity of the voting model were 0.8755, 0.9239, 0.9192, and 0.6923, respectively (Fig. 3A), and the overall performance of the model exceeded that of previously reported models^17,18,40 (Supporting Information Table S14). In addition to the 50% cell viability threshold, 80% cell viability prediction models were developed to identify compounds with weak toxicity. The RF model performed the best, with an AUROC of 0.8408 and accuracy, precision, recall, and specificity metrics of 0.8112, 0.8225, 0.8910, and 0.7875, respectively (Fig. 3B). Identification of highly safe substances from massive samples in the early stage of drug discovery could significantly reduce subsequent costs in further development, which can be accessed either by choosing nontoxic samples or by excluding toxic substances. Among the metrics, “precision” indicates the proportion of the predicted nontoxic substances that are truly nontoxic, i.e., how likely the chosen nontoxic substance is. The “specificity” represents the ratio of truly predicted toxic samples to total toxic samples, i.e., the ability of the model to exclude toxic substances from all samples. Both optimal models, the voting model for 50% cell viability prediction and the RF model for 80% cell viability prediction, performed well in cell viability prediction, with precision scores of 0.9239 and 0.8225, and specificity scores of 0.6923 and 0.7875, respectively. Besides classical machine learning algorithms, deep learning algorithms are capable of extracting features directly from raw data to enable end-to-end prediction and offer flexible architectures for stable and scalable training with large datasets. Therefore, four DL algorithms, including multilayer perceptron, convolutional neural network TabResnet, and TabNet, were also utilized to develop prediction models and the AUROCs of the four models ranged from 0.4851 to 0.8751 (Supporting Information Table S15). The unfavorable performance of DL models might be due to the high-dimensional data (over 5000 features) and relatively small size of samples (1243 compounds), which made the models overfit.

The generalizability of the optimal models across cell lines was evaluated in the HEK293 (kidney), HepG2 (liver), and HT-29 (colon) cell lines. Gene expression profile data for the cells treated with the test compounds and relevant cell viability data (Supporting Information Tables S16–S18 and Fig. S3) were used for this evaluation. Both the voting model (50% cell viability) and the RF model (80% cell viability) exhibited good performance (Supporting Information Table S19), with AUROC values ranging from 0.7898 to 0.9318 and 0.7622 to 0.8146, respectively, across the three cell lines (Fig. 4), indicating the two optimal models could be used for cell viability prediction across cell lines. In addition to AUROCs, “precision” and “recall” are also important in this case as they represent the ability to select nontoxic substances, with “precision” representing the proportion of the predicted nontoxic substances that are truly nontoxic, and “recall” representing the ratio of truly predicted nontoxic samples to the total number of actual nontoxic samples. As shown in Table S19, the precision and recall scores of two optimal models ranged from 0.6969 to 0.9217 and from 0.8023 to 0.9833, respectively, indicating that the models have a strong ability to predict nontoxic substances and could contribute to the identification of “safe” substances for drug development. Additionally, regression models also were developed, and their Pearson correlation coefficients ranged from 0.7025 to 0.7422 (Supporting Information Fig. S4). In summary, these findings indicate that, by using the models built in this study, cytotoxicity at different thresholds can be accurately predicted in diverse cell lines, which could aid in primary cytotoxicity screening across cell lines from massive amounts of transcriptome data.

Interpretability is pivotal because it provides information on the decision-making process of the model and uncovers elemental features. To determine the importance of the Feature Genes for cytotoxicity classification, two optimal models for 50% and 80% cell viability prediction, the voting model and the RF model were interpreted using SHAP analysis. SHAP is a feature attribution method that estimates the contribution of each feature to the prediction. The absolute value of SHAP represents the extent of feature attribution, while positive and negative SHAP values indicate the contributions of the feature to the prediction of noncytotoxicity and cytotoxicity, respectively. In this study, the SHAP values among samples in the training set were calculated (Supporting Information Tables S20 and S21), and the average absolute SHAP values across samples were used to measure the global importance of features in the model (Fig. 5A and B). Additionally, the distributions of the SHAP values of the Ensemble Features among the samples in the training set were displayed, which revealed the contribution of each Feature Gene to a certain sample (Fig. 5C and D). As expected, for both the voting model (50% cell viability prediction) and the RF model (80% cell viability prediction), the top 10 features with the highest feature importance were involved in cell cycle events (CCNA2, CDK1, and RUVBL1 in the voting model; TOP2A, CCNA2, CDK1, POLE2, MCM10, TMPO, and CENPX in the RF model; Fig. 5A and B). Of the remaining differences between the two models due to cytotoxicity level variance, the top 10 genes of the voting model (50% cell viability prediction) were involved in metabolic and immune pathways (ATP8B4, DERA), and the top 10 genes of the RF model (80% cell viability prediction) were involved in the cellular stress response (TSC22D3, HES1) (Fig. 5A and B). In addition, the distributions of individual gene expression patterns and SHAP values among the samples were analyzed. As shown in Fig. 5E and F, the expression level of the top 10 features of the two optimal models, the voting model for 50% cell viability prediction and the RF model for 80% cell viability prediction was correlated with the direction of cytotoxicity/viability. The results showed that the interpretation of the prediction models provided insights into cytotoxicity, including the quantification of feature importance and the direction of regulation. As important features made vital contributions to models, two re-fit models (the voting model for 50% cell viability prediction and RF for 80% cell viability prediction) were developed with their top 10 features, and the AUROCs of the two re-fit models across four cell lines are shown in Fig. 6. Compared with the AUROCs of the models with all ensemble features, the re-fit models also exhibited good performance (0.77 to 0.87 by re-fit voting vs. 0.79 to 0.93 by voting model; 0.79 to 0.84 by re-fit RF vs. 0.76 to 0.81 by RF model), indicating that the important features obtained by SHAP analysis had substantial impacts on cytotoxicity/viability prediction. Notably, the good performance and generalization of the 10-features models across cell lines may provide an opportunity for cytotoxicity/safety prediction with minimal cellular transcription information, which will make the models widely applicable.

As high-dimensional data, cellular transcriptomic data represent molecular changes that embody valuable information for compound pharmacological mechanism studies. However, for NTI compounds, the median therapeutic concentration is close to or even overlaps with the cytotoxic concentration, resulting in the co-occurrence of biofunctional activities and toxic effects that obscure the primary mechanisms involved. To exclude confounding effects, the CTS genes, which represent the major gene group causing cytotoxicity, were measured and defined. First, a total of 41 genes, which had higher mean absolute SHAP values than the average absolute SHAP value of 228 Feature Genes in the RF model, were acquired as the most relevant subset of genes associated with cytotoxicity (Fig. 7A; Supporting Information Table S22). The variable genes correlated with each gene in the subset were subsequently measured with a correlation coefficient greater than 0.5 as the threshold, and 545 genes (Fig. 7A; Supporting Information Table S23) were defined as CTS genes. Subsequently, to investigate the application of CTS genes in the MOA study, two known NTI compounds, cyclosporine A and triptolide, were analyzed by using transcriptome data with or without CTS gene removal to verify the contribution of these CTS genes to the MOA study of NTI compounds.