The overall workflow of our study is shown in Figure 1. Initially, the overlap of drugs, included in each of the gene expression cell data sets, was investigated using VennDetail (Guo and McGregor, 2020) to create a Venn pie chart showing the various drug testing subsets across the six cell lines (Figure 2). Each of the non-gene expression datasets (FAERS, MOLD2, and TOX21) were treated as individual datasets, while the gene expression data were merged across cell lines to build classifier models. In general, we used standard preprocessing techniques, including removing zero variance features and missing values. DILI1 and DILI3 suffered from class imbalance (Table 3). For all non-gene expression data, to mitigate this issue, we attempted three oversampling techniques, including synthetic minority oversampling technique (SMOTE) (Chawla et al., 2002), random oversampling examples (ROSE) (Menardi and Torelli, 2014), and a random upsampling of the minority classes. SMOTE balances data by randomly creating artificial samples between two nearest-neighbor samples, while ROSE uses a smoothed bootstrap technique to resampled the data (Chawla et al., 2002; Menardi and Torelli, 2014). For comparison, models were built using imbalanced data as well. Before training non-gene expression datasets, they were standardized to have a mean of zero and a standard deviation of one. Preprocessing details specific to each dataset as well as some characteristics of the data are discussed below.

The CAMDA organizers provided us with FAERS data for all 617 drug compounds. Of these, 422 were grouped as “training data”. This dataset contains 20 features corresponding to information on the percentage of reported adverse events for each drug compound by gender and age group demographic. After removing highly correlated features, we upsampled the data to cater to the class imbalance by randomly sampling with replacement from the minority class to balance the majority class. An additional preprocessing step was to create two new variables, namely “male ratio” and “female ratio”, taking into account all reported events irrespective of the gender, all reported DILI events irrespective of gender, and the percentage of reported DILI events by gender.

In addition to the FAERS dataset, we were provided with concentration-response information of 600 drugs. Of these, 412 were designated as “training data”. Thirty-two features corresponded to concentration-response curve ranks. Out of all 412 drugs for training, 57 drugs were removed for missing values. In addition, we removed highly correlated features using an arbitrary cutoff of 0.82 and catered to the class imbalance by using SMOTE.

Alongside the FAERS and TOX21 data provided, we had access to the 2D molecular descriptors or structural information of these 617 drug compounds. 422 of these drugs were designated for training. There were 777 features for each drug compound with each feature corresponding to MOLD2 descriptors. To cater to class imbalance, we upsampled minority classes, as well as ROSE, and SMOTE.

The L1000 assay data used in this study is a high-throughput gene expression assay that measures mRNA transcript abundance of 978 landmark genes based on an inference algorithm to infer the expression of 11,450 additional genes in the transcriptome (Subramanian et al., 2017). Utilizing simulation, it has been observed that this reduced representation of the transcriptome can recapitulate around 80% of the relationships of measuring the entire transcriptome directly. In this study, 12,328 deidentified predictors genes were provided by the CAMDA organizers with Z scores to indicate transcript abundance. The treatment time and dosage of each drug were selected by the CAMDA committee to produce the largest available dataset for both test and training data.

Since not all drugs were tested in each cell line data made available, we utilized the Kruskal-Borda (Kru-Bor) merging algorithm in the GeneExpressionSignature R package (Li et al., 2013). This approach allowed us to generate a unified drug-induced expression signature across cell types since many drugs were not tested in the PHH or HepG2 liver cell lines. The Kruskal algorithm (Kruskal, 1956) finds a minimum spanning forest of an undirected edge-weighted graph while the Borda merging method (Saari, 2000b; 2000a) uses ranked options in order of preference to determine the outcome. Thus each closest neighbor in rank merges one by one until a unified signature is formed. Following merging, the top and bottom 100, 250, 500, and 1,000 ranked genes were selected as drug signatures for feature selection.

A method of feature selection utilized across the merged signatures produced via our Kru-Bor merging was based on a gene’s significance (p-value < 0.01) in predicting the DILI score via a Fisher’s exact test. If a gene is included in the top or bottom 100, 250, 500, or 1,000 ranked list, depending on the model data, for any drug it would be assigned a one (True), or if it fell outside of that range it would be assigned a zero (False). The classifier for each type of DILI was also one (DILI positive) or zero (DILI negative). We used these classifiers to identify if these highly influenced genes were predictive of a drug being DILI positive or DILI negative with a p-value cutoff of 0.01.

The prediction of DILI was treated as a binary classification problem for each DILI type. That is, for each of DILI1, DILI3, DILI5, and DILI6, outcomes were split between “positive” and “negative”. We used a 5-folds cross validation repeated 100 times, and a random search strategy to search for the best parameters for each model. The data was made available such that training and test sets had been pre-identified. Importantly, we did not have access to the correct labels for the test data. Models were built using traditional machine learning algorithms within the caret (Kuhn, 2008) package in R version 4.0.0 (R Core Team (2020), 2020).

The machine learning algorithms we used are suitable for classification tasks. They include a Logistic Regression (LR) (Cox, 1958), Linear Discriminant Analysis (LDA) (Li and Jain, 2009), Decision Trees (DT) (Quinlan, 1986), Support Vector Machines (SVM) (Cortes and Vapnik, 1995), Naïve Bayes (NB) (Hand and Yu, 2001), a One-layer Neural Network (Nnet), and a Random Forest (RF) algorithm. LR and LDA are generally categorized as linear classification models, with an assumption that the data follows a normal distribution. Given a set of predictors, LR aims to build a linear model of these predictors by minimizing the sum of squared residuals. LDA uses the prior probability of belonging to a class to estimate posterior probabilities by using Bayes’ Theorem. DT and RF are often classified as trees and rules-based algorithms. Given a set of predictors, a decision tree works by using if-else conditions to build a definitive set of rules using splits. The challenge usually lies in determining optimal situations to apply a “then”-clause (or a split). In RF, similar conditional statements are used. However, instead of using the entire sample of data for tree-building, RF uses many independent subsamples from the training data to build small decision trees. Each small decision tree classifies an observation by voting. Neural networks and SVMs are generally grouped as non-linear algorithms. Neural networks (in our case, a multilayer perceptron i.e. a neural network with one hidden layer), are modeled after how neurons in the human brain work. The outcome or prediction is a linear combination of the hidden layer(s) transformed by a non-linear activation function. There are several activation functions used, depending on whether the problem is a regression or classification problem. In our case, we used a sigmoidal or logistic function, since we were dealing with a classification problem. SVM aims to find support vectors or data points that separate the different classes as much as possible. Intuitively, these data points are the most difficult to separate (the reasoning is that they lie very close to one another and to the hyperplane or decision boundary), and are thought of to be important in separating classes. There are different flavors of SVMs depending on the kernel used (kernels are similar to non-linear activation functions used in neural networks). In the current study, we used polynomial, linear, and a radial basis function kernels.

To evaluate the performance of our models, we focused on the area under the ROC curve (AUC) value as well as the specificity, sensitivity, accuracy, and MCC of the models on the test set. The AUC value is a widely-used metric in binary classification problems. An AUC value of one indicates a perfect classifier, i.e. a model that is perfectly able to separate both classes, while an AUC value of 0.5 indicates a model that predicts at random. Depending on the application domain, AUC values of 0.7 and above are usually acceptable. Specificity measures the ratio of negative classes that were correctly identified by the model out of all negative classes, while sensitivity measures the ratio of positive classes that were correctly identified by the model out of all positive classes. These metrics are affected by how the target labels are structured and passed to the algorithm, and they range from 0 to 1. Additionally, we evaluated the performance of our models on the test set by calculating the balanced accuracy of prediction. Balanced accuracy is the average of the sensitivity and specificity or the average of the fraction of correct labels that are predicted correctly (by the model) within each class. We used this metric because we observed that there was class imbalance within our datasets regardless of DILI type.

MCC is particularly useful in datasets of different class distributions (or imbalanced data) because it considers all of the false and true positives and negatives. It is calculated from the confusion matrix of a model and its values range from +1 to −1, with +1 indicating a perfect classification, 0 indicating random classifications, and −1 indicating no relationships between the observed and predicted classes.

In an attempt to improve the classification accuracy of our models, we used three ensemble voting approaches, namely soft voting, hard voting, and a weighted voting approach. These ensemble methods work best when there are varying algorithms of different strengths i.e. algorithms having varying underlying assumptions about the data, and when each one has reasonable predictive power (Kuncheva, 2002; Van Erp et al., 2002). Using the gene expression data provided by CAMDA 2018 organizers, Sumsion and others (Sumsion et al., 2020) used hard and soft voting ensemble methods in an attempt to improve prediction accuracy on DILI risk. As an extension of their work, we hypothesized that since we have access to larger and more diverse datasets, we could capture different aspects of predicting DILI types and use these ensemble methods to improve prediction.

Hard voting, also known as majority voting, takes into account the predicted class labels of each classifier (or voter) (Ruta and Gabrys, 2005). Voting is done by counting how many class labels (for each class) were predicted among all classes. The class label with the highest count is taken to be the predicted class label for that observation. On the other hand, soft voting considers the probabilities of each class label by each classifier (Lin et al., 2003). In other words, it considers how certain each classifier is about the class labels. For each class label, the probabilities are averaged, and the label with the highest average probability is taken as the predicted class label for the observation.

The third approach to voting involves using a weight to skew predictions towards the most certain models. In our approach, we used the AUC of each classifier as a weighting parameter for the output probabilities. This was done to take into account that some classifiers might have better predictive power and should be given preference in determining the outcome of the voting. To weigh each probability, we multiplied the probabilities of each predicted class by the AUC and divided this by one subtracted from the weight, that is, the AUC of that model. Afterward, weighted probabilities were treated just as in soft voting: by taking the average of all resulting weighted probabilities belonging to each class. The class label with the higher average was taken as the predicted class for that observation. Therefore, the predicted class, y ^ , of observation, given an output set of class membership probabilities across many models, P , is given by: C ( y ^ | P ) = a r g m a x     1 m ∑ i = 1,2 … m m ( w i m ∗ p i c w i m − 1 )   Where m is a model, w i m is the weighting parameter for a model, p i c is the probability of class membership.

The performance of FAERS data in predicting each of the DILI types can be seen in the bar plots in Figure 3. While we built many models, we compared and picked the best three models based on the AUC values to predict DILI class on the test set. We noticed that using the raw data (without resampling), models achieved classification accuracy between 0.51 and 0.55 and MCC between 0.04 and 0.14 on the training set and did not do noticeably better on the test set (accuracy: 0.49 to 0.59, MCC: −0.03–0.22). On the other hand, using resampled datasets improved the accuracy of the models on the training set to a range of 0.61–0.94 (MCC: 0.47–0.89). Using these models to predict the DILI class of the test set showed a slight improvement in the accuracy (0.52–0.62). The MCC, however, was between 0.04 and 0.24.

The top three models built using TOX21 data (using the AUC as the criterion) were evaluated on the test set (Figure 4). Using the data as is, without resampling, the accuracy of the training data was between 0.50 and 0.57 (MCC: −0.02–0.17). As expected, the models failed to generalize to the test set (accuracy: 0.50 to 0.59, MCC: −0.04–0.19). Again, we observed that resampling slightly improved the accuracies of these models on the training set (accuracy: 0.62–0.76, MCC: 0.25–0.54). Yet, there was no major improvement on the test set (accuracy: 0.50 to 0.58, MCC: −0.01–0.20).

Similarly to how the FAERS data was handled, we selected the top three performing models built using MOLD2 data in each category (resampled or non-resampled) to predict the DILI class of the test data (Figure 5). Models built using the non-resampled MOLD2 dataset gave accuracies of 0.50–0.54, showing that the models were randomly predicting the classes (MCC: 0.00–0.17). This performance was similar on the test set (accuracy: 0.50 to 0.66, MCC: −.01–0.36) with a slight improvement. Similarly to what we observed using FAERS data, resampling the dataset improved both the accuracy and the MCC of the training set (accuracy: 0.71–0.78, MCC: 0.56–0.76) but could not generalize better than non-resampled MOLD2 data to the test set (accuracy: 0.51 to 0.67, MCC: 0.14–0.36).

Cellular RNA expression levels in the form of microarray data have been previously investigated for their ability to predict DILI with limited predictive power (Chierici et al., 2020). In the current study, the L1000 data from the Connectivity Map was used including both the measured landmark genes as well as the inferred transcriptome. We built models using each expression data to investigate which cell lines were most successful in predicting DILI. Table 4 summarizes the model results based on our training data. However, due to the limitation of each cell only providing expression response data from a subset of drugs (Figure 2) involved in the training and test data, accuracy based on test data was not meaningful. Additional processing steps for this data involved merging across the six cell lines to generate a representative signature, testing different cutoffs for the amount of highest- and lowest-ranked genes to utilize, as well as a feature selection for determining predictor genes.

The models built using the merged expression signatures with the highest AUC from the training data were evaluated on the test set. The training and test results are summarized in the bar plots in Figure 6. None of the cell expression signatures performed well when predicting DILI3, DILI5, or DILI6 with an accuracy ranging from 0.39 to 0.64 and MCC values ranging from −0.03 to 0.1. These models did have some limited success predicting DILI1 with the merged SVM 1000 model performing the best, reaching an accuracy of 0.67 but an MCC of 0.10 (Supplementary Table 1). The poor predictability of DILI3 status by these models was unexpected with the accuracy of the best model being 0.49 with 0.33 sensitivity and 0.66 specificity. The limited success in predicting DILI5 and DILI6 was expected based on the positive and negative control construction of these DILI classes, which are not reflected in the gene expression data.

Since the top three individual models did not perform well on the test set (Supplementary Table 1), we asked if aggregating the top three models in an ensemble approach could improve the accuracy. To test this, we applied three ensemble voting methods namely soft voting, hard voting, and weighted voting. Hard voting gave accuracies of 0.39 and 0.37 on DILI1 and DILI3, respectively, while soft voting gave an accuracy of 0.44 and 0.40 for DILI1 andDILI3, respectively (Figure 7). Soft voting slightly improved the accuracy of these models most likely because it considers membership probabilities rather than predicted class labels. We observed that weighted voting slightly improved the accuracy: 0.54 for DILI1 and 0.60 for DILI3. Our weighted approach considers both the probabilities and the AUC of the models and emphasizes the contribution of models with higher AUCs. Sumsion and others used similar approaches (soft and hard voting) with gene expression data resulting in decreased accuracies (Sumsion et al., 2020). Compared to their study, our approach improved the accuracies of the models. However, our method(s) does not report MCCs because we do not have access to the true positives, true negatives, false positives, and false negatives in the test data.