Anticancer drug synergy prediction based on CatBoost

Datasets

We used the datasets from reference (Fan, Cheng & Li, 2021), which considered the NCI-ALMANAC (https://wiki.nci.nih.gov/display/NCIDTPdata/NCI-ALMANAC) dataset as the source of anticancer drug synergy data. NCI-ALMANAC contains data of 59 cancer cell lines and we only considered drugs that have at least one target gene (68 drugs). A total of 130,182 samples were used for model training and testing. All of the NCI-60 cancer cell line characteristics (expression, mutations, copy number variant (CNV), etc.) could be obtained from CellMinerCDB (https://discover.nci.nih.gov/rsconnect/cellminercdb/) (Luna et al., 2021). Drug target data was provided by DrugBank (Wishart et al., 2018), and the drug molecular properties could be processed using the RDKit package in Python. Cell line characteristics generally include gene expression profiles, mutation numbers, copy numbers, etc., while drug characteristics usually include morgan fingerprints, drug target information, monotherapy information, etc. In this work, we used gene expression profiles as cell line features, and morgan fingerprints, drug target information, and monotherapy information as drug features. The shape of the dataset is 130,182 samples with 2,064 columns of features, of which 470 are cell line features and 1,594 are drug features.

We used Python’s pandas package to process the data. Firstly, pandas is used to extract four columns of data (drug1, drug2, cell line, score) required in this article, which are drug1, drug2, cell line and synergy score of anticancer drugs. This is a matrix with 130,182 rows and four columns, which is named as the master matrix in this article. Next, through the constructed master matrix, the cell line features are constructed in this article. The gene expression profile contains the gene expression profiles of 470 genes under 59 cell lines, from which a 59 × 470 matrix was constructed in this article.

Then, a zero matrix with 130,182 rows and 470 columns was constructed. Then, the cell lines in the first row of the main matrix were used to get the corresponding row of the gene expression profile matrix, and then assigned to the first row of the zero matrix. The construction of drug features is roughly the same as above, and since each group of samples contains two drugs, the drug features of the two groups of drugs need to be merged. After obtaining all the cell line features and drug features, the NumPy concatenate function is used to concatenate the data by column (axis = 1) to obtain the complete feature data matrix.

CatBoost

Based on the constructed features, this work used a new GBDT framework named CatBoost to predict drug synergy (Fig. 1). CatBoost, as an oblivious trees-based learner, has fewer parameters, and supports for categorical variables, as well as high accuracy GBDT framework. Compared with other methods, CatBoost does not require manual one-hot encoding when processing categorical data, which greatly simplifies the data preprocessing process. In addition, CatBoost reduces the risk of overfitting by using a technique called “ordered boosting” and improves the generalization ability of the model by using a symmetric tree structure. Existing methods, such as XGBoost and LightGBM, although also very powerful and popular, may require additional preprocessing steps when dealing with categorical features, such as manual feature encoding. Moreover, these algorithms may require more hyperparameter tuning to achieve the best performance. As for why CatBoost was chosen over LightGBM, although LightGBM may outperform CatBoost in training speed, CatBoost performs better in terms of generalization accuracy and area under the curve (AUC) on a specific dataset, although this difference may be small. CatBoost’s native support for categorical features and its ability to reduce model overfitting make it potentially a better choice in some cases. In addition, the design of the algorithm structure of CatBoost helps to reduce the gradient boosting bias and prediction offset, thus improving the accuracy of the model. Algorithm 1 gives the pseudo-code of the GBDT algorithm.

Compared with the traditional decision trees, the oblivious trees can better deal with the continuous attribute and the classification problem in the high dimensional space, which improves the prediction ability of the model. Since the oblivious trees can completely ignore some attribute information, it is very small for the existence of noise and missing values in the data. In addition, the oblivious trees use the principle of symmetry, and through the restriction of depth and number of branches, it can effectively avoid the problem of overfitting, and improve the model’s generalization ability. CatBoost also uses the concept of ordered boosting, a permutation-driven approach, which trains the model on a subset of the data while computing the residuals on another subset, thus it can prevent target leakage and over fitting. Algorithm 2 gives the pseudo-code of the ordered boosting algorithm.

In a nutshell, CatBoost has the following peculiarities: 1. excellent performance, it surpasses most advanced machine learning algorithms in terms of performance; 2. robustness, it reduces the need for much hyperparameter tuning and decreases the chance of overfitting, which makes the model more versatile; 3. practicality, it can handle both categorical and numerical features; 4. scalability, it supports customized loss functions.

The code and datasets of this work could be obtained from: https://github.com/AnnaGuan/CatBoost (DOI: 10.5281/zenodo.14189227).

Experimental setup

To make CatBoost generalizable to unseen datasets, we used stratified k-fold cross-validation experiments to test the model. Stratified k-fold cross-validation is an enhanced version of k-fold cross-validation for unbalanced datasets. In stratified k-fold cross-validation, the whole dataset is divided into k equal-sized copies, and the positive and negative ratio of the label variable in each fold is the same as the percentage in the whole dataset. This work used stratified 5-fold cross-validation to evaluate the model performance. The dataset of all samples was divided into five equal and unique parts, of which four parts were treated as training data and the remaining one as test data. Each part was regarded as test data in turn, and its synergy scores were calculated using the learned model.

By using the stratified 5-fold cross-validation, we also adjusted the parameters of CatBoost to obtain the most optimal model. The main parameters we adjusted included the iterations in the range of [400, 500, 600, 700, 800, 900, 1,000], the depth in the range of [5, 6, 7, 8, 9, 10, 11], and the learning rate in [0.05, 0.5, 0.1]. By comparing the model prediction results under each parameter combination, the optimal parameter combination to train the model was set as: iterations = 600, depth = 9, learning rate = 0.1.

Drug synergy prediction usually has two ways: one is regression task and another is classification task. In this article, we have considered both of these two pattens. When conducting synergy prediction as a classification task, the synergy scores were binarized, i.e., synergy labeled as 1 and antagonism as 0. Since there are many samples with synergy scores close to 0, it is crucial to choose an appropriate threshold to binarize the synergy score. In this work, we optimized the threshold by implementing stratified 5-fold cross-validation to achieve optimal balance, and the value of threshold was finally set as 10.

To evaluate the prediction performance of CatBoost, we compare its synergy prediction ability with the following models: 1. deep neural network (DNN), a deep learning method based on feedforward neural networks; 2. XGBoost, a classical integrated learning method based on gradient boosting; 3. logistic regression, a log-linear regression model based on the sigmoid function. The hyperparameters used in these models were determined by stratified 5-fold cross-validation.

The following statistical measures and methods are mainly used in this article: 1. receiver operating characteristic area under the curve (ROC AUC): it is mainly used for classification tasks to measure the performance of a classifier. 2. Precision area under the recall curve (PR AUC): this is used to evaluate the performance of a classification model, especially when the dataset is imbalanced. 3. mean squared error (MSE): a measure of the average of the squared differences between the predicted and actual values, used for regression tasks. 4. Pearson correlation coefficient (PCC): this is used to assess the linear relationship between the predicted and actual values.

All the experiments were done on Windows 11 on a work computer with the CPU model: 11th Gen Intel (R) Core (TM) i5-11300H @ 3.10 GHz and GPU: Intel (R) lris (R) Xe Graphics.

Performance evaluation

We measured the prediction performance of our model in terms of the receiver operating characteristic area under the curve (ROC AUC) and the area under the precision-recall area under the curve (PR AUC) for classification. To further characterize the prediction performance of CatBoost, we also provided typical performance metrics for regression tasks: mean square error (MSE) and Pearson correlation coefficient (PCC). After carrying out stratified 5-fold cross-validation based on the optimal parameter setting mentioned above, we can compute the values of four metrics for our model and other methods for comparison. The results were shown in Table 1, from which we found that CatBoost outperformed the other three models on all metrics. For example, by comparing the ROC AUC values, CatBoost achieves 0.9217, improved 0.0099 over the next best model DNN; by comparing MSE, CatBoost obtained 0.1365 that was improved 6% over the next best model XGBoost.

To further investigate the quality of CatBoost predictions, we compared the distribution of synergy scores on cell lines between the predicted result and the actual values (seen in Fig. 2). Panel (a) shows the box plots of the distribution of the actual synergy scores for each cell line, and panel (b) shows box plots for the synergy scores predicted by CatBoost. Each point represents the actual or predicted synergy score for a drug combination. To understand how well our model predicted synergy for drug combinations on individual cell lines, we also gave the average PCC between predicted and actual synergy scores for each cell line (Fig. 2C). The order of the cell lines in Fig. 2C is the same as in panel (a) and (b). From panel (c), we could observe that the correlation between all cell lines ranged from 0.52 to 0.83.

Additionally, by calculating the PCC, we analyzed the performance of CatBoost from two viewpoints, anticancer drug and cancer cell line, respectively. Figure 3 gives the PCC distribution over 68 drugs, where each bar represents the average PCC of all combinations of the specific drug tagged below the bar. The color of the bars corresponds to the drug target. As shown in Fig. 3, the PCC values for all anticancer drugs ranged from 0.51 to 0.88, where 31% of the drugs have PCC values below 0.6% and 39% above 0.7. We also found that the color of the bars is scattered, which means there is no clear association between the PCC and the drug target. Therefore, the differences of performance between drugs could not be explained by target-based mechanisms. Moreover, the number of drugs acting on the same target does not affect the performance.

The PCC distribution over cell lines is illustrated in Fig. 4. Each bar is the average PCC between actual and predicted synergy scores for a cell line. The cell line names are shown on the x-axis. The color of the bars corresponds to the type of tissue from which the cell line originated. From this figure we can see that PCC values for all cell lines vary from 0.52 to 0.83, with only three cell lines having values below 0.6. Among them, more than 44% of the cell lines have PCC above 0.7. Furthermore, the colors of the bars are scattered and no correlation between tissue type and PCC can be observed. Therefore, the performance differences could not be explained by the tissue type of cell lines either.

Based on the above analysis, this article argues that the difference in prediction performance is not clearly related to the tissue-specific mechanism and the tissue type/target coverage in the dataset. These differences may be due to the fact that some biological mechanisms are better modeled than others. This may be related to the measurements obtained from cell lines in this article and whether they are able to capture these biological processes.

Feature interpretability

While it is important for a model to predict drug combinations with synergistic effects, it is also important to explain why the model could predict synergistic drugs. Explainable machine learning has slowly become an important research direction in machine learning in the past few years. SHAP is a model interpretation package developed in Python, which can interpret the output of machine learning models. Inspired by cooperative game theory, SHAP constructs an additive explanatory model in which all features are considered “contributors”. For each prediction sample, the model generates a predicted value, and the SHAP value is the number assigned to each feature in the sample, which can be calculated as follows:

(1) $$SHAPvalues=\sum _{s\subseteq M\left\{i\right\}}{\displaystyle \frac{\left|S\right|!\left(\left|M\right|-\left|S\right|-1\right)!}{\left|M\right|!}\times f\left(S\cup \left\{i\right\}-f\left(S\right)\right)}$$where $M$ is the set of features, $S$ is a subset of $M$ excluding feature $i$, and $f\left(S\right)$ is the prediction of the model on the feature set $S$.

For binary classification, we can verify the SHAP calculation using the following formula:

(2) $$basevalues+\sum SHAPvalues\left[0\right]=\ln \left({\displaystyle \frac{p}{1-p}}\right).$$

In this work, the SHAP values of all sample features were calculated by SHAP, and the 100 features with highest scores were obtained (see Table S1). We found that among the top 100 features, monotherapy information ranked first, suggesting that monotherapy information was more useful for predicting anticancer drug combination synergy. Statistically, 88 out of these 100 features belong to drug, of which Morgan’s molecular fingerprint takes the largest proportion. The remaining 12 important features were genes in gene expression profiles, including PTK2, CCND1, GNA11, CRKL, ERBB2, WNT2B, CTBP2 etc. In order to investigate the importance of cell line features and drug features, we implemented a comparative experiment through using only cell line features or only drug features respectively for model training and predicting, the results of which are shown in Table 2. We can see that the results of the model using only cell line features are much worse than those using only drug features, both of which are not as good as the original results. Apparently, drug features play a more prominent role in drug synergy prediction than cell line features. The conclusion we reached is consistent with the findings of Janizek, Celik & Lee (2018).

For the genes ranking by the front, we have made a case study to further investigate the significance of our model. For instance, PTK2, also known as FAK, was found to associate with tumors and functions primarily in the inhibition of apoptosis (Yi et al., 2021). The upregulation of FAK expression, including high protein expression as well as overactivation, is present in almost all tumor tissues, such as lung, gastric, colorectal, uterine, and melanoma cancers. According to studies (Chuang et al., 2022; Xie et al., 2022), patients with high FAK expression have a significantly shorter survival than those with lower FAK expression, making FAK a promising indicator for early diagnosis of cancer and for predicting the prognosis of patients.

The gene CCND1, or cell cycle protein, plays a role in regulating CDK kinase (Akervall et al., 2003). This protein interacts with the tumor suppressor protein Rb that could positively regulate the expression of this gene. Mutations, amplification and overexpression of this gene alter the course of the cell cycle, which are frequently observed in a variety of tumors, and may contribute to tumorigenesis (Hosokawa & Arnold, 1998).

ERBB2, or erb-b2, is a 185-KDA cell membrane receptor encoded by the proto-oncogene erbB-2 (Penuel et al., 2001). Clinically, erb-B2 expression is strongly correlated with patient prognosis, and patients with high erb-B2 expression are prone to tumor metastasis and have short survival. Because of the significant difference in expression levels between normal and tumor cells, erb-B2 has become an ideal target for tumor immunobiotherapy and is currently a hot molecule in the field of tumor therapy research (Mcdonagh et al., 2012).

To sum up, the top-ranked genes all play important roles in cancer, and such roles are also extended to the field of anticancer drug combination prediction and captured by the model in this article, making an important contribution in the prediction. Obviously, the prediction results of the model in this article have certain biological significance.

Feature reduction

In order to further verify the impact of the front-rank features with high SHAP values on the model efficacy, we selected the top 400 features as input features and performed a new prediction using CatBoost, which was written as new CatBoost. The new CatBoost experiment has the same settings as the previous CatBoost experiment, with the only difference being the number of features used. Considering the balance between computing speed and performance, we chose the top 400 features instead of less, since too few features would result in performance degradation. This feature selection method can be applied to other data sets. The verification experiments were conducted both in stratified 5-fold cross validation and stratified 10-fold cross validation. The prediction performance of new CatBoost is given in Table 3. By comparison, it is clear that when only these 400 features are taken into account, the performance of CatBoost is not affected much. Compared with the original CatBoost, the results of new CatBoost are even slightly improved and the results of the stratified 10-fold cross validation are better than that of the stratified 5-fold cross validation.