Prediction of cancer cell line-specific synergistic drug combinations based on multi-omics data

Framework of XDDC

The process of XDDC is illustrated in Fig. 1. We employed drug combination-cell line data from the NCI-ALMANAC database and collected relevant drug and cell line feature data from various public databases. To mitigate the sparsity of the input matrix, we conducted data filtering and feature compilation. Subsequently, we utilized five machine learning methods and a deep learning model for modeling, fine-tuned the models and performed evaluations. Finally, predictions were made using the optimal model, and the results were validated.

Drug and pathway data screening

We collected the names of 104 drugs from NCI, obtained the target information of drugs from HCDT (Chen et al., 2022), and screened 340 human-related pathways (including genetic data on pathways) from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa & Goto, 2000). Considering that the sparsity of the matrix will affect the predicted results, we first screened the drug and pathway data according to the following criteria: (1) the drug has at least 10 target genes; (2) there are at least 10 genes in the pathway; and (3) at least half of the drugs in each pathway can act on it. As there was little data about the direct action of drugs on pathway, we assumed that the drug act on a pathway when the target gene of the drug is in the pathway. As a result, we retained 87 drugs and 157 pathways.

Drug features

For drugs, we collected information based on the chemical structure, toxicity, and pathways associated with drug targets. Chemical fingerprinting is the most commonly used structural feature for drugs, allowing the representation of chemical properties at the molecular level (Carracedo-Reboredo et al., 2021). We utilized fingerprints to represent drugs as binary vectors indicating the presence (1) or absence (0) of various chemical features (such as S-S, N-O, and C = CN bonds). We collected the SMILES of drugs from HCDT and PubChem (Kim et al., 2021), generated 166-length MACCS molecular descriptors, and represent each drug as a binary 0–1 vector using the RDkit library (Lovric, Molero & Kern, 2019) in Python. For each drug combination, we used the sum of the two drug vectors as a feature. (i.e., 0 indicates neither drug possesses the structural feature, one indicates only one drug in the combination has the feature, and two indicates both drugs in the combination have the feature). We retained only 109 descriptors as the structural features of drugs after ensuring that each MACCS molecular descriptor was present in at least ten drugs.

We collected toxicity features from the SIDER database, with 16 drugs having no matching results. For these 16 drugs, we gathered information from DrugBank (Wishart et al., 2018) based on an overview of drug toxicity. Keywords were extracted from the descriptions of toxicity in DrugBank and matched to the characteristics in SIDER, capturing toxicities for these drugs. We selected toxicity labeled as LLTs in SIDER and ensured that each side effect had at least one corresponding drug, resulting in the retention of 97 features of drug toxicity. Similarly, each drug was represented as a binary vector, and the feature for each drug combination was represented as the sum of the two drug vectors.

Out of the 87 drugs selected from NCI, a total of 1,571 target genes were obtained from the HCDT database. We filtered out 157 human biological pathways in KEGG, generating 157 features for the subsequent analysis. The individual drug-pathway matrix was denoted as D.

(1) $$\begin{array}{c}D=\left[\begin{array}{c}{d}_{11}\\ ⋮\\ d_{m1}\end{array}\begin{array}{c}{d}_{12}\\ ⋮\\ d_{m2}\end{array}\begin{array}{c}\dots \\ \ddots \\ \dots \end{array}\begin{array}{c}{d}_{1,n}\\ ⋮\\ d_{mn}\end{array}\right],d_{i,j}=\{\begin{array}{c}0\\ 1\end{array}\begin{array}{c}\mathrm{\#}\left(d_i\cap p_j\right)=0\\ \mathrm{\#}\left(d_i\cap p_j\right)>0\end{array}\mathrm{\#}.\end{array}$$

The drug-pathway matrix D represented whether a drug acted on a particular pathway, where $d_{i,j}$ was computed based on the intersection of the target gene set $d_i$ of drug i and the gene set $p_j$of pathway j. A value of 0 indicated that the drug did not act on the pathway, while a value of 1 indicated that the drug acted on the pathway.

Due to the interactions between pathways, we utilized the Jaccard coefficient to assess the correlation between pathways, and we denoted the pathway-pathway matrix as P:

(2) $$\begin{array}{c}P=\left[\begin{array}{c}{p}_{11}\\ ⋮\\ p_{n1}\end{array}\begin{array}{c}{p}_{12}\\ ⋮\\ p_{n2}\end{array}\begin{array}{c}\dots \\ \ddots \\ \dots \end{array}\begin{array}{c}{p}_{1,n}\\ ⋮\\ p_{nn}\end{array}\right],p_{i,j}={\displaystyle \frac{\mathrm{\#}\left(p_i\cap p_j\right)}{\mathrm{\#}\left(p_i\cup p_j\right)}\mathrm{\#}}.\end{array}$$

The P matrix was an n*n matrix, where $p_{i,j}$ was calculated based on the gene sets of pathway i and pathway j. $\mathrm{\#}\left(p_i\cap p_j\right)$ indicated the number of genes in the intersection of pathway i and j, $\mathrm{\#}\left(p_i\cup p_j\right)$ indicated the number of genes in the union of gene i and j.

Considering that a drug could act on multiple pathways and that there were interactions between pathways, we refined the effect score of a drug on a pathway according to Formula (3):

(3) $$\begin{array}{c}\tilde{D}=DP=\left[\begin{array}{c}{\tilde{d}}_{11}\\ ⋮\\ {\tilde{d}}_{m1}\end{array}\begin{array}{c}{\tilde{d}}_{12}\\ ⋮\\ {\tilde{d}}_{m2}\end{array}\begin{array}{c}\dots \\ \ddots \\ \dots \end{array}\begin{array}{c}{\tilde{d}}_{1n}\\ ⋮\\ {\tilde{d}}_{mn}\end{array}\right],{\tilde{d}}_{ij}={\overrightarrow{D}}_i^T\cdot {\overrightarrow{P}}_j=\sum _{k=1}^n{d}_{ik}\cdot p_{kj}\mathrm{\#}.\end{array}$$

The improved drug-pathway matrix $\tilde{D}$ was calculated by multiplying matrices D and P, where ${\overrightarrow{D}}_i^T$ represented the ith row vector of matrix D and ${\overrightarrow{P}}_j$ represented the jth column vector of matrix P. In addition to assessing the effects of drugs on pathways, the $\tilde{D}$ matrix also considered the interactions between pathways based on their similarities.

The final drug combination-pathway matrix A was obtained from the drug-pathway matrix D and $\tilde{D}$:

(4) $$\begin{array}{c}A=\left[\begin{array}{c}{a}_{1,1,1}\\ ⋮\\ a_{m,m,1}\end{array}\begin{array}{c}{a}_{1,1,2}\\ ⋮\\ a_{m,m,2}\end{array}\begin{array}{c}\dots \\ \ddots \\ \dots \end{array}\begin{array}{c}{a}_{1,1,n}\\ ⋮\\ a_{m,m,n}\end{array}\right],a_{i,j,k}={\displaystyle \frac{1}{2}\left(d_{i,k}\cdot {\tilde{d}}_{j,k}+d_{j,k}\cdot {\tilde{d}}_{i,k}\right)\mathrm{\#}}.\end{array}$$where $a_{i,j,k}$ represented the impact of the combined use of drug i and drug j on pathway k. In the end, the matrix A obtained through three-dimensional matrix operations not only considered the effects of drug combinations on pathways but also considered the mutual interactions between pathways.

Multi-omics data processing for cell lines

The multi-omics data describing cancer cell lines, including EXP (expression), MUT (mutations), CNV (copy number variations), METHY (methylation), and RNAi (RNAi-related data), were downloaded from the Cancer Cell Line Encyclopedia (CCLE) (Nusinow et al., 2020) and the Dependency Map (DepMap 22Q2) portal (http://depmap.org).

The level of gene expression reflected the level of gene activity in a cell and was a key link between functional genes and phenotypic changes. Gene expression data was obtained by RNA sequencing and processed using the RSEM tool to obtain Transcripts Per Million (TPM) expression values. We obtained 19,221 genes across 1,406 cell lines after a log2(TPM+1) transformation.

Gene mutations led to the loss or gain of protein function, directly affecting cell function and the occurrence and development of diseases, helping to identify key disease-causing genes and potential targets for targeted therapies. The mutation data was derived from the collective output of MuTect1, MuTect2, and Strelka, encompassing information on all somatic point mutations and indels called in DepMap. This dataset contained mutation information for 18,784 genes in 1,771 cell lines.

Copy number variation affected gene expression and function by altering gene dose, which might lead to tumor formation or the development of other complex diseases. Copy number variation scores were calculated by mapping genes to the segment-level calls and computing a weighted average along the genomic coordinate. There were copy number variation information for 25,368 genes in 1,766 cell lines after log2(CN ratio + 1) transformation.

DNA methylation referred to the addition of methyl groups to DNA molecules and was usually associated with the inhibition of gene expression. It played an important role in cell differentiation, development and disease. Methylation data was obtained by quantifying CpG islands using reduced representation bisulfite sequencing. Following Wang et al. (2020), we first removed genes that showed NA in more than 70% of cell lines. Subsequently, we utilized the K-nearest neighbor (KNN) to impute the remaining NA values. Finally, we calculated the average methylation value for each gene across all cell lines, resulting in data for 20,183 genes and 842 cell lines.

RNA interference mechanisms regulated gene expression by targeting specific mRNAs and was used to identify key genes that drove changes in cell function. RNAi data processing followed a similar approach to methylation data. Genes with NA values in more than 70% of the cell lines were first removed. Subsequently, the KNN was applied to impute the remaining NA values. The final dataset contained RNAi information for 12,500 genes and 710 cell lines after calculating the mean value.

There were 38 common cell lines between DepMap and NCI ALMANAC.

Cell line-pathway features

To provide a more comprehensive interpretation of the biological mechanisms of the model, we used pathways as a bridge connecting drug features with cell line features, which required the construction of a cell line-pathway matrix. For the four sets of cell line features (expression, copy number variation, methylation, and RNAi), we utilized the Gene Set Variation Analysis (GSVA) (Hanzelmann, Castelo & Guinney, 2013) method to obtain cell line scores for pathways. We employed the Jaccard coefficient to determine the overlap between genes with mutations and the gene set of pathways for each cell line. This process resulted in a mutation-related cell line-pathway feature matrix. Each set of cell line omics data results in 157 features.

NCI-ALMANAC data preprocessing

The NCI-ALMANAC dataset comprises 5,356 drug combinations, 104 drugs, 60 cell lines, and nine types of cancer. We first excluded all relevant data about single-drug actions. Then, since both NSC753082 and NSC761431 represent the drug vinorelbine, we removed the data related to NSC753082.

The drug combinations in NCI-ALMANAC were treated at different concentrations (3*3 or 5*3) in the same cell line (Xia et al., 2018). We averaged the synergy scores at different concentrations for each drug combination and cell line, named as CombScore.

(5) $$\begin{array}{c}Combscore={\displaystyle \frac{1}{n}\sum _{a,b}{E}_{a,b}-P_{a,b}\mathrm{\#}}.\end{array}$$

We did not consider the concentration of drugs A and B. The mean difference between the expected growth score and the observed percent growth score was used as the CombScore.

We defined the drug pairs in each cell line as synergistic or not based on the CombScore. We categorized the CombScore using a threshold of 4 and compared it to the previous study using a threshold of 10 (Li et al., 2020). Drug pairs with CombScore greater than the threshold were considered synergistic, otherwise they were considered non-synergistic.

Drug combination modeling with machine learning and deep learning algorithms

We utilized five machine learning models—logistic regression, lasso regression, random forest, support vector machine (SVM) and eXtreme Gradient Boosting (XGBoost), along with one deep learning model, TabNet (Joseph, Joseph & Prasad, 2022). These models were trained on the pre-processed feature sets to evaluate their predictive performance on drug combination synergy.

Model adjustment and validation assessment

To evaluate the predictive performance of all machine learning models, we randomly divided the entire dataset into training and test sets with an 80:20 ratio. By using the StandardScaler from the scikit-learn library for feature normalization, we ensured that the features of both the training and test sets were adjusted to follow a standard normal distribution with a mean of 0 and a standard deviation of 1.We employed the GridSearchCV function from scikit-learn to perform a grid search and select the optimal parameters for improved performance. To address the issue of an imbalanced sample distribution, we employed techniques such as adjusting sample prediction weights, undersampling, and oversampling during the model development process, aiming to rebalance the ratio of negative to positive samples to 2:1. The sampled dataset was further divided into five equal-sized subsets for five-fold cross-validation to avoid problems such as oversampling after data processing and to ensure the stability of model predictions. Evaluation metrics including accuracy, precision, recall, F1-score, area under the receiver operating characteristic curve (AUC), and area under the precision-recall curve (AUPR) were used to evalute the binary classification models. Due to the imbalance between negative and positive samples, we focused primarily on AUPR values to assess the predictive ability of the models, emphasizing the correct prediction of true positive samples.

Drug validation experiments

We cultured T-47D cells (human mammary ductal carcinoma cells) purchased from Pricella and used them within 2 weeks of receipt. The cells were characterized by short tandem repeat (STR) profiling. The cell lines were passaged until growth reached 70–85% of the flask. The passaging process involved washing the residual culture medium with phosphate-buffered saline (PBS), digesting the cells with tryptic cell digest, and terminating the digestion with T-47D cell-specific medium (Pricella, CM-0228), which was prepared as a cell suspension. Passaging was performed every 48 h at a ratio of 1:3. During cell processing, fresh culture medium was prepared and cell counting was performed by diluting the cell suspension to 50,000 cells/ml. The single-cell suspension was then spread and 200 μL of single-cell suspension was added per well in a 96-well plate, resulting in 10,000 cells per well. The cells were incubated in a cell culture incubator at 37 °C, 5% CO₂, and 95% humidity for 24 h (Pettersen et al., 2007).

For drug treatments, all drugs were purchased from Selleck. For single drug studies, cells were treated with six concentrations of bortezomib (0.5, 2.5, 10, 20, 30, and 40 nmol/L), with each concentration repeated eight times. In addition, cells were treated with four concentrations of bleomycin (50, 100, 200, and 400 nmol/L), and the experiment was also repeated eight times for each drug concentration. Cell survival was determined after 48 h and the optimal concentration was selected based on the rate of inhibition of cell growth by each drug. To study drug combinations, eight replicate treatments were performed with two drugs added to the cell lines at the optimal concentration, and the experiment was repeated twice to account for the order of drug administration. The results of the two experiments were recorded as Combination 1 and Combination 2.

Finally, cell viability assays were performed using the CCK-8 method by aspirating the drug-containing medium and adding 100 μL of fresh medium (Cai et al., 2019). To each well, 10 μL of CCK8 activity assay reagent was added and incubated for 1 h at 37 °C, 5% CO₂ and 95% humidity. Cell activity was detected using a microplate reader at 450 nm. The cell activity formula was calculated as follows:

(6) $$\begin{array}{c}Cellviability\left(\mathrm{\%}\right)=\left[A_d/A_n\right]\times 100\mathrm{\%}\mathrm{\#}.\end{array}$$

We compared the effectiveness of single drugs and drug combinations using the method of Li et al. (2020):

(7) $$\begin{array}{c}IDAcombscore=\mathrm{△}via-\mathrm{△}via\times H{R}_{C/Mbest}\mathrm{\#}\end{array}$$where $\mathrm{△}via$ was equal to the mean viability when cell lines were treated with the best monotherapy minus the mean viability when cell lines were treated with the drug combination. $\mathrm{H}{\mathrm{R}}_{\mathrm{C}/\mathrm{M}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$ was the hazard ratio between the drug combination efficacy and the best monotherapy efficacy.

Exploration of NCI-ALMANAC dataset

The NCI-ALMANAC dataset was one of the largest known drug combination datasets, containing 104 drugs, 5,356 drug pairs and 60 cell lines across nine cancer types, with each cancer type containing 2–9 cell lines. To mitigate the impact of data sparsity on the input features, using a Combscore threshold of four as an example, we selected 87 drugs with at least 10 targets each, resulting in 3,741 drug pairs, some of which have been experimentally validated by the NCI (see Figs. 2A, 2B). When the threshold is set to 10, specific details of drug combinations are illustrated in Fig. S1. We selected 38 cell lines with comprehensive information on expression, mutation, copy number variation, methylation, and RNAi (Figs. 2C, 2D). In total, 128,445 drug-drug-cell line data points were retained for model training, with an additional 13,713 data points reserved for testing (see Table S1). Simultaneously, we downloaded 340 human-relevant pathways from KEGG, filtering for those in which at least half of the drugs acted and where each pathway contained at least 10 genes. This resulted in a final set of 157 pathways.

Due to the limited number of positive samples, the sparsity of the input feature matrix may affect the accuracy of predicting synergistic effects. Therefore, we applied specific criteria to filter drugs, pathways, and cell lines (see Methods for details). The feature representation of each sample is a 1,148-dimensional vector. After normalizing the input features, we employed various models for training, including linear regression, lasso regression, random forest, support vector machine, XGBoost and TabNet, and comparing their results.

Model selection and evaluation

To compare the model performance under two different thresholds, we used all feature sets (1,148 features) as inputs for model training. With a threshold of 4, there are 14,371 positive samples (drug-drug-cell line combinations with synergistic effects) and 114,075 negative samples. With a threshold of 10, there are 3,575 positive samples and 124,871 negative samples. Due to the imbalanced distribution of the classification results in the samples, we adopted weight assignment to increase the probability of predicting positive samples, with the default being ‘class_weight = balance’.

Comparing the results of the six models under the two thresholds (Table 1), it could be observed that the models with a threshold of 4 had higher precision, recall and AUPR values (Fig. 3A), whereas the models with a threshold of 10 had higher accuracy and AUC values (Fig. 3B). Considering the imbalance in the distribution of positive and negative samples, we finally chose the threshold with higher AUPR values, i.e., a threshold of 4, as the criterion for determining whether drug combinations show synergistic effects.

By comparing the various evaluation metrics of the six models, the overall performance of the XGBoost model outperformed the other models. The performance and training time of XGBoost is significantly influenced by the choice of parameters. To improve the accuracy of the model predictions, we used grid search to determine the values of several commonly used variables in XGBoost: ‘learning_rate’ was set to 0.1, ‘max_depth’ was set to 7, and ‘n_estimators’ was set to 700. In addition, we utilized a five-fold cross-validation method to prevent over- or under-fitting of the models on certain datasets, with the best AUC and AUPR being 0.865 and 0.564 respectively. We continued to use lattice search and cross-validation in subsequent model adjustments to improve the predictive performance of the models.

Feature selection

We tested our model on different combinations of feature categories to select the optimal input feature set. There were three types of drug features, including MACCS molecular descriptors, drug toxicity, and pathways, and five types of cell line features, including expression, copy number variation, mutation, RNAi, and methylation. We randomly combined these eight feature sets into a total of 255 input feature sets. Subsequently, we utilized XGBoost with the same optimal parameters for modeling and evaluated the results using the AUPR value (details in Table S2).

The evaluation results, displaying the top 5 AUPR values and the remaining combinations of multiple feature sets, were presented in Table 2. Considering that recall, F1-score, and AUPR were commonly used metrics to evaluate imbalanced data classification, we ranked these three metrics separately and calculated their average ranks. The table shows that the model using a combination of seven feature sets—MACCS, toxicity, pathway, expression, CNV, RNAi, and methylation—performed the best, achieving an AUPR value of 0.5604 and ranking first. It can also be observed from Figs. 4A–4C that it demonstrated a high F1 score and Recall values. On the other hand, the model using all eight feature sets ranked only eighth in terms of AUPR.

Meanwhile, we observed that using a single feature set, whether it’s drug related or cell line related, the predictive performance is in the bottom 50% across all models. Specifically, when using only cell line features, the precision, recall, and F1 values all equate to 0. This may be attributed to the small sample size of 38 cell lines compared to the total samples of 14,371 positives and 114,075 negatives, causing the predictions to heavily favor negative samples, indicating no synergistic effect of drug combinations on cell lines. However, when combining drug and cell-line features, the predictive performance of most models improves. Taking drug pathway features as an example, the inclusion of methylation data increases the ranking of the model by 21 positions compared to using cell line pathways alone. In addition, using multiple cell line data is more effective than adding a single cell line data. For example, modeling with only drug pathways and methylation data performs less effectively than models incorporating drug pathways along with CNV, methylation, and mutation data. Furthermore, we observed that using only drug features is less effective than combining drug and cell line feature sets. Taking drug toxicity and drug pathways as an example, their predictive performance is inferior to the feature set that includes CNV data. In addition, introducing more cell line features slightly improves the predictive performance of the model.

In summary, we ultimately chose the feature set consisting of MACCS, drug toxicity, drug pathways, gene expression, copy number variation, RNAi, and methylation as the input features for the model. With this feature set, we achieved the best evaluation results using XGBoost.

Correction of unbalanced data

To address the issue of imbalanced sample distribution, we employed the following three methods to re-model the previously selected optimal feature set combination. (R1) Adjusting positive and negative sample weights: Given the larger proportion of negative samples, we adjusted the weights based on the ratio of negative to positive sample quantities. This adjustment aids convergence in the context of imbalanced samples. (R2) Undersampling: Using the imblearn package from sklearn, we reduced the number of positive samples to achieve a balanced ratio of 1:2 between positive and negative samples. (R3) Oversampling: Employing the SMOTE method, we synthesized samples for the minority class, maintaining a balanced ratio of 1:2 between positive and negative samples. The evaluation results of establishing models using the three methods are shown in Table 3.

Among the three methods, the oversampling approach consistently outperformed the weighting and undersampling approaches in terms of various evaluation metrics. Although the method of weight allocation can yield a better ACC value, but its prediction results tended to favor negative samples, which was contrary to our goal of obtaining more positive samples. By referring to Figs. 5A, 5B, in the ROC and PR curves of AUC and AUPR, the areas under the curves for the XGBoost model after oversampling processing reached 0.9685 and 0.9595, respectively, demonstrating excellent predictive performance. The ROC and PR curves for the other two methods were shown in Fig. S2.

In addition, as our focus was more on the predictive performance of positive samples, we aimed to achieve better recall, F1-score and AUPR during the model evaluation process. Through a horizontal comparison of multiple models under the same method, the results (see Figs. 5C–5F) indicated that regardless of the method used to handle imbalanced data, the XGBoost model exceled in the classification task of predicting synergistic effects. Compared to other classical machine learning models commonly used for such tasks, such as logistic regression, lasso, SVM, and random forest, the XGBoost model showed a significant improvement in various evaluation metrics, including notable enhancements in accuracy and precision (see Figs. S3E–S3F). Therefore, we ultimately decided to adopt the XGBoost model trained with oversampling as the preferred training model.

Feature importance evaluation

To better explain the biological mechanisms behind drug synergies, we connected information on drug combination effects to cell line information via biological pathways. We selected five feature sets related to pathways from the seven feature sets and calculated the importance of each feature (details in Table S3). These included one drug feature (drug pair pathway) and four cell line features (expression-pathway, methylation-pathway, CNV-pathway, and RNAi-pathway). From the A, it could be observed that the contributions of these five feature sets were relatively similar, accounting for 47% of the total feature importance. Despite the drug MACCS fingerprint being the most important feature set, acknowledged as an excellent feature for studying drug synergies, it lacked relevant support for understanding the biological mechanisms of drug synergy in our research.

To investigate the most important pathway information among the five feature sets, we performed a more detailed analysis. We selected the top one, two, three, four, five, and 10 features from each feature set, repeated the modeling using the XGBoost model, and plotted the trends of various evaluation metrics (Fig. 6B). We found that using the top 10 features in each feature set, ranked by contribution, as input data can essentially achieve the same results as training with all features.

Summarizing the top ten features (pathways) from each feature set (Fig. 6C), we observed that the pathways associated with each feature set were almost unique. Only the linoleic acid metabolism pathway ranked in the top ten for both CNV and RNAi features, and the prion disease pathway was prominent in METHY and RNAi. To understand the functional information associated with these pathways, we annotated them based on the primary functions in KEGG. Most of the top ten pathways for each feature set were related to human diseases functions. However, there were still distinctions in emphasis among feature sets. For instance, DDP still includes pathways related to environmental information processing, whereas CNV features were associated with organismal systems. In particular, among the top three ranked features (Table S4), DDP pathways were mainly related to environmental information processing, EXP was predominantly linked to organismal systems, CNV, and METHY primarily involved pathways related to human diseases, and RNAi showed a broader range of pathways. This also indicated that different omics information provided distinct biological pathway insights, enhancing the mechanistic scope of biological pathways through multi-omics information.

Prediction and validation in unlabeled data

In NCI-ALMANAC, only partial experimental data on the interaction of drug combinations with cell lines are provided. We processed the remaining drug combination information in the same way, including 113,713 drug combination-cell line data as a validation set. We used the trained model to make predictions, and the predicted results were shown in Table S5. Among these, there were 287 positive results, indicating drug combinations predicted to have a synergistic effect on cell lines. Figures 7A–7D illustrated all synergistic drug combination-cell line pairs predicted for breast cancer, leukemia, non-small cell lung cancer, and ovarian cancer, respectively. Results for other cancers are shown in the Fig. S3.

To validate the authenticity of the predicted results, we conducted a literature search on all positive outcomes and selected one for experimental verification. Three cases were validated through the literature. Specifically, the combination of doxorubicin hydrochloride and celecoxib showed a synergistic effect on the HL-60 (TB) cell line. The combination of 4′-epiadriamycin and erlotinib hydrochloride inhibits the growth of the A549 cell line, and Arsenic trioxide in combination with cisplatin acts on the ovarian cancer-related IGROV1 cell line.

He et al. (2015) research indicated that in HL-60 cell lines, both MDR1 and COX2 mRNA were both expressed at low levels. After treatment with 0.02 μg/ml daunorubicin for 72 h, the expression of MDR1 and COX2 mRNA increased by 3.68 and 3.29 times, respectively. Subsequent treatment with celecoxib revealed that the induction of MDR1 mRNA expression by daunorubicin was significantly upregulated, while COX2 mRNA showed a slight decrease. Overall, compared to daunorubicin alone, the combination with celecoxib resulted in a 3.07-fold decrease in the effective inhibitory concentration (IC₅₀) for HL-60 compared to daunorubicin alone, and the use of celecoxib alone did not inhibit the growth of HL-60 cell lines.

Similar to the combination of daunorubicin and celecoxib, He et al. (2016) suggested that the combination of daunorubicin and erlotinib can inhibit the growth of A549 cell lines associated with non-small cell lung cancer. Additionally, research by Sweatman & Israel (1987) indicates that adriamycin (ADR) and 4′-epiadriamycin (epi-ADR) showed comparable biological properties in most in vitro and in vivo test systems. Moreover, 4′-epiadriamycin has lower cardiac toxicity. Therefore, we infer that the combination of 4′-epiadriamycin and erlotinib can also inhibit the growth of A549 cell lines. Additionally, Zhang et al. (2009) found that compared to using cisplatin or arsenic trioxide alone, the combination of cisplatin and arsenic trioxide can inhibit the growth of ovarian cancer-associated IGROV1 cell lines at lower doses. The dose reduction index (DRI) ranges from 1.75 to 5.21 under different concentrations.

Drug validation experiments

At the same time, we carried out drug experiments on the synergistic inhibitory effect of bleomycin and bortezomib on breast cancer cell line Tmur47D. We measured the cell proliferation activity of bleomycin at four concentrations and bortezomib at six concentrations for 48 h and determined that the optimal concentration of bleomycin was 200 nmol, and that of bortezomib was 40 nmol/L. The results of the drug experiments in Table 4 showed the effect of the combination of two drugs on cell activity at the optimal concentration through multiple repeated experiments. The combination experiments for the two drugs were conducted repeatedly, employing the studied drug combination efficacy calculation method of Ling & Huang (2020). Using the cell activity detected for each single drug and the cell activity of the combination drug, the obtained IDAcomboscore values were 0.014 and 0.018, respectively, both exceeding 0.004 (see Fig. 7E). Therefore, these two drugs were considered to have synergistic effects.