2.1 Data

Drug screening data and the transcriptomic data of the cell lines were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) [11] and Cancer Cell Line Encyclopedia (CCLE) [12]. The molar concentration of a drug needed for half inhibition of cell growth (IC50) as well as molecular information of cell lines in GDSC and CCLE, such as gene expression profiles, mutation, and copy number variation were downloaded using PharmacoGx R package [13]. The max concentration values for drugs were obtained from GDSC website http://cancerrxgene.org.

GDSC contains 439 drugs and 1124 cell lines; however, some IC50 values were not inserted in GDSC. Therefore, we purified the data using some pre-processing steps similar to several previous works [5, 9, 14, 15]. After applying data pre-processing steps, we obtained 555 cell lines and 98 drugs.

Moreover, CCLE dataset contains the information for 24 drugs. It should be noted that we require the max concentration values for drugs. Therefore, according to the previous studies [4], we restrict the CCLE dataset on the drugs that their max concentration values can be obtained from GDSC website. The final CCLE dataset, after applying pre-processing steps, contains 363 cell lines and 18 drugs.

In addition to the mentioned information, the PaDel descriptors of drugs in PubChem [16] were extracted using the rdkit package in Python [17]. The target proteins for drugs were gained from GDSC and DrugBank [18] databases as well as literature. The interaction information for proteins and drugs are obtained from STRING [19] and STiTCH [20] databases.

2.2 Pre-processing

We construct three types of matrices for assessing the gathered data. The first matrix is the response matrix which contains the drug response information. The rows of the response matrix corresponds to the cell lines and its columns are related to the drugs. The elements of response matrix are either the IC50 values or sensitivity labels. The second type of matrices is cell line feature matrices, in which the rows pertain to the cell lines and the columns represents different features. The last type of matrices is drug feature matrices, in which each row represents a drug and drug features are organized in the columns. After constructing these matrices, pre-processing steps were applied on them in order to impute the missing values, remove samples with significantly low information, calculate the cell line similarities and drug similarities, and make the data suitable for the proposed method. The pre-processing procedure includes the following four steps:

Imputing missing values

Converting IC50 values into binary categories

Similarity calculation

Standardization and normalization

These steps are elaborated in the following subsections.

2.2.1 Imputing the missing values

There are numerous missing values in the IC50 response matrix, cell line feature matrix based on copy number variation, and cell line feature matrix based on mutation profile. In order to remove the samples that have a significant lack of information, we omitted the drugs that have missing IC50 for more than half of the cell lines. Moreover, we removed the features of cell lines that have the missing values for the majority of cell lines. Afterward, we excluded the cell lines with missing entries for more than half of the columns in each of these matrices. Nevertheless, there were still some missing values in the matrices; therefore, we imputed these missing entries using a weighted mean of other entries. If we did not omit the samples with a significant lack of information and tried to impute all missing entries first, the obtained information would not be much reliable. The percentage of imputed data must be low enough to maintain the authenticity of data. For example, in the case of GDSC dataset, the raw IC50 matrix contains about 49% missing values. After removing the drugs and cell lines with a great extent of missing values, the obtained IC50 matrix contains only 2.7% missing pairs. Imputing such a limited fraction of data does not damage the data authenticity.

We used the following strategy to impute the remaining missing values, which is similar to the imputation procedure used in previous studies. [5, 14]. Let E(c_i) be the gene expression profile of cell line c_i and I(c_i, d_j) be the IC50 value for cell lines c_i against drug d_j. The missing value for I(c, d) is imputed as the Eq 1.

where ||.||₂ is the norm-2 and N(c) is the set of indices for cell lines which are the k nearest neighbors of the cell line c with respect to the distance function D. We considered k = 10 for imputing missing values in GDSC and CCLE.

Moreover, the missing entries in copy number variation and mutation matrices are imputed similarly. Let V(c, g) and M(c, g) be the copy number variation and mutation status of gene g in cell line c. The missing values for V(c, g) and M(c, g) are imputed according to Eqs 3 and 4, respectively.

It is noteworthy that gene expression profiles of cell lines are considered for calculating sample distance because the cell line feature matrix based on gene expression profiles does not contain any missing value. Therefore, the distance function D can be computed for all pairs of cell lines.

2.2.3 Similarity calculation

Previous studies have frequently confirmed that similar cell lines yield similar responses to similar drugs [5–7, 25]. Therefore, the machine learning approaches can learn to predict the drug response using the similarities between cell lines and the similarities between drugs. Three types of cell line similarities and three types of drug similarities were computed for GDSC and CCLE datasets. The similarity calculation procedure is elaborated in the following.

The genomic features of cell lines are mainly characterized using gene expression profiles, copy number variation, and mutation profiles. The gene expression similarity between c_i and c_j cell lines is represented by SC_E(c_i, c_j), which is computed by Pearson Correlation Coefficient (PCC) between the gene expression profiles of the c_i and c_j (See Eq 5).

where E(c_i, g) denotes the expression of gene g in cell line c_i and $$ \overline{E}\left(c_i\right)$$ represents the mean expressions of all genes in cell line c_i.

Moreover, SC_V represent the cell line similarity based on copy number variation. Let V(c_i, g) be the copy number variation of gene g in cell line c_i and $$ \overline{V}\left(c_i\right)$$ mean of copy number variations of all genes for cell line c_i. The similarity of cell lines corresponding to copy number variation is calculated as Eq 6.

In addition to the mentioned similarities between cell lines, another important feature in cancer cell lines is mutation profiles. Gene mutations play crucial functions roles in cancer development and progression [26]. Suppose M(c_i, g) be a binary value, showing the mutation status of gene g in cell line c_i, where it equals “1” if gene g is mutated in cell line c_i and “0” if it is wild type. SC_M(c_i, d_j) represent the mutation similarity of cell lines c_i and c_j which is defined based on Jaccard Index (JI) of their mutation profiles.

One of the frequently used similarities for drugs is the similarity of chemical substructures because chemical substructure of a drug determines its functionality up to a good extent [27–29]. The chemical substructure similarity of two drugs d_i and d_j is shown by SD_S(d_i, d_j) which is calculated using the JI of PaDel descriptors of drugs d_i and d_j. Let P(d_i, l) be the lth element of the PaDel descriptor for drug d_i. The SD_S(d_i, d_j) value is computed as Eq 8

Another informative similarity of drugs can be obtained based on the interaction of drugs with other chemicals and proteins. The STiTCH database provides a comprehensive resource which presents the relationships between chemicals and proteins [20]. The relations in STiTCH are based on various sources such as experimental evidence from ChEMBL [30] PDSP Ki database [31], Protein Data Bank (PDB) [32], pathway databases including KEGG [33], Reactome [34], and NCI nature pathway interaction database [35] as well as other databases such as DrugBank [18] and MATADOR [36]. In addition to the experimental evidences from reliable databases, it uses text mining, experimentally biochemical data, gene fusion, and genomic context prediction [20]. Hence, the network obtained from STiTCH database provides an extensive insight about the drugs. The STiTCH network for CCLE drugs is illustrated in Fig 1. The oval nodes in this figure represent the drugs in CCLE and other adjacent drugs, while the circle nodes indicate the neighbor proteins. The STiTCH network for GDSC drugs is presented in the S2 File (S1 Fig in S2 File). The weights of edges between nodes in STiTCH network for GDSC and CCLE drugs are presented in S1 and S2 Tables in S1 File, respectively.

Fig 1

STiTCH network for drugs in CCLE.

We computed drug similarities based on STiTCH network (SD_N) according to the Eq 9, where N_i represents the immediate neighbors (both protein neighbors and chemical neighbors) of drug d_i. Therefore the JI of the neighbor nodes in STiTCH network is used as the second similarity of drugs.

In addition to the above drug similarities, we calculated another similarity based on PPI network of target proteins. To this aim, we obtained drug targets and considered the target proteins for all drugs in the dataset as the target protein set. The target protein set for GDSC and CCLE are presented in S3 and S4 Tables in S1 File, respectively. Afterwards, we acquired PPI network for the target protein set from STRING database [19]. Fig 2 represents the STRING PPI network for target protein set of CCLE drugs. The STRING PPI network for GDSC drugs is presented in the S2 File (S2 Fig in S2 File). It is noteworthy that STRING database uses various types of data such as Gene ontology terms, pathways experimental evidence, and text mining features to compute the interactions between proteins. Hence, the weights of edges in STRING PPI network are the combination of various evidence. The weights of edges between nodes in STiTCH network for GDSC and CCLE drugs are presented in S5 and S6 Tables in S1 File, respectively.

Fig 2

STRING network for drugs in CCLE.

To compute the PPI-based similarity of drugs (SD_P) we constructed a bipartite graph for each pair of drugs. To explain clearly, suppose that TP_i and TP_j are the sets of target proteins for drugs d_i and d_j, respectively. Bipartite PPI graph G(i, j)^PPI is constructed such that the set of nodes in one part is TP_i, the set of nodes in another part is TP_j and the edges are defined based on STRING PPI network. For example, assume that that TP_i = {P_i1, P_i2, P_i3} is the target proteins of drug d_i, while TP_j = {P_j1, P_j2, P_j3, P_j4} is the target proteins of drug d_j. The bipartite graph G(i, j)^PPI is illustrated in Fig 3. The weights PPI edges between two parts in this graph are set by the weights of PPI in STRING network. Afterwards, we applied the maximum weighted matching algorithm [37] on this bipartite graph and consider the summation of weights of matching edges as the PPI-based similarity of drug pairs. This similarity shows the extent of accordance between the set of target proteins of two drugs. Therefore, SD(d_i, d_j) is a high value, if the set of target proteins for d_i, d_j match highly.

Fig 3

Schematic representation of the bipartite graph between the set of target proteins for a drug pair.

2.2.4 Standardization and normalization of similarity matrices

After computing all similarities, we standardized the similarity matrices in order to ensure that all similarity values range from 0 to 1. Since all similarities were computed based on PCC or JI, the similarity values range from -1 to 1; therefore, this standardization transforms the values in the range of [0, 1]. To clearly explain the standardization process, let S(i, j) be an entry in a similarity matrix. Its standardized value is represented by $$ \widehat{S}(i,j)$$ and is computed according to Eq 10.

It is notable that performing standardization on similarity matrices does not change the sorting of distances between samples because this is a linear transform.

In the next step, we normalized the standardized similarity matrices using the symmetric normalized Laplacian [38]. The symmetric normalized Laplacian is a well-defined transform of the similarity matrix with several favorable algebraic and spectral characteristics such as being positive definite and diagonally dominant [39]. Moreover, its prevalence use in other problems such as spectral clustering [40] and drug target interaction prediction [41] justifies that using symmetric normalized Laplacian matrix represents the similarity of samples and shows the structural properties in a better way [42]. For each similarity matrix S, the normalized similarity matrix $$ \mathcal{S}$$ is obtained as Eq 11.

In this equation, D is a diagonal matrix and D_i,i = ∑_j S_i,j. All diagonal elements in D are non-zero. Hence, D^−1/2 is a diagonal matrix and its entries are reverse values of the square root of elements in D. It should be noted that this normalization does not affect the correlation between various types of similarity matrices. On the other hand, applying this normalization improves the model speed and leads to quicker convergence of the CDSML.

2.3 Applying classification method

The classification process is inspired from ADRML method [9], but with binary response matrix B, rather than IC50 value matrix I. Moreover, it uses a threshold to convert obtained results into binary labels. This method has several steps, including:

Decompose B_m×n matrix into X_m×k and Y_n×k, such that B ≈ XY^T

(a)Initialize X⁽⁰⁾ and Y⁽⁰⁾ randomly.

(b)Compute a loss function defined as the summation of mean square error (MSE), similarity conservation terms and regularization terms.

(c)Update X and Y matrices according to the Newton’s method.

(d)Repeat steps 1.b and 1.c until X and Y matrices converge.

Decompose B^T into two latent matrices W and Z similar to step 1.

Compute $$ \tilde{B}={\displaystyle \frac{1}{2}}(X{Y}^T+W{Z}^T)$$ as the predicted sensitivity matrix.

Use a threshold to convert $$ \tilde{B}$$ into a binary matrix $$ \widehat{B}$$.

These steps are shown in Fig 4 and explained elaborately in the following.

Fig 4

The overall framework of the proposed method.

The features cell lines, IC50 and C_max values were obtained from GDSC. Besides, drug substructures were downloaded from PubChem, STiTH network from STiTCH datbase, and PPI network from STRING database. The sensitivity associations between cell line-drug pairs were specified and used as the objective of manifold learning. The cell line similarities and drug similarities were calculated, standardized, and normalized. These similarities were considered as the regularization terms in manifold learning. The output of manifold learning is a predicted score matrix that assigns cell line-drug pairs into sensitive resistant classes by a threshold.

2.3.1 Compute loss function

To classify the cell line-drug pairs, B_m×n is decomposed into two latent matrices X_m×k (cell line latent matrix) and Y_n×k (drug latent matrix) with lower rank. The decomposition must satisfy the following constraints:

Decomposition Mean Square Error (MSE): The multiplication of latent matrices XY^T must be an appropriate estimation of binary response matrix B.

Regularization: The elements of X and Y matrices should not grow excessively.

Cell line similarity conservation: The similar cell lines must have not too far latent row vectors in X matrix.

Drug similarity conservation: The similar drugs must have not too far latent row vectors in Y matrix.

Considering all the above constraints, the loss function is defined according to Eq 12.

where α and β are the regularization and similarity conservation coefficients. X(i) and Y(i) denotes the ith rows in X and Y matrices, respectively. The symbol SC in Eq 12 can be substituted by SC_E, SC_C, or SC_M. Two latent matrices X and Y were updated using Newton’s method to minimize the loss function iteratively. X⁽⁰⁾ and Y⁽⁰⁾ were initialized randomly and afterwards, X^(k) and Y^(k) were updated using the rules defined in Eqs 13 and 14, respectively.

These matrices were updated until they do not change significantly in two subsequent iterations. Specifically, when ||X^(k+1) − X^(k)|| + ||Y^(k+1) − Y^(k)|| < 0.01, the convergence criteria is met. The detailed formulae for updating latent matrices are described in S2 File.

2.3.2 Predicting sensitivity/resistant labels

As described above, we decomposed B into two latent matrices and after convergence, the last estimated latent matrices X^(k) and Y^(k) were multiplied to estimate B:

Then, B^T is decomposed into W_m×k and Z_n×k using the described method. This is done due to the fact that the predicted labels for samples (c_i, d_j) and (d_j, c_i) must be equal. After the convergence, B^T is estimated using the following equation:

The predicted sensitivity matrix ($$ \tilde{B}$$) was calculated as the average of $$ {\tilde{B}}_1$$ and $$ {\tilde{B}}_2$$. It should be noted that the estimated matrix $$ \tilde{B}$$ is not binary valued; thus, we converted it to a binary-valued matrix $$ \widehat{B}$$ using a threshold. So that the cell line-drug pairs were assigned to sensitive/resistant classes. The computation of best threshold is explained in Section 2.4.

2.4 Evaluation criteria

The predictive performance of models was assessed using stratified five-fold cross-validation on cell line-drug pairs. To this aim, the set of all cell line-drug pairs were partitioned randomly to five subsets of equal sizes such that the fraction of sensitive over resistant pairs was almost equal in all subsets. Four subsets were considered as the training data and the evaluation criteria were computed on the remaining subset. This procedure was iterated for each subset and the criteria were averaged over the iterations. The stratified five-fold cross-validation was repeated 30 times in order to prevent bias in partitioning the dataset. The most prevalent evaluation criteria for classification problems are defined in the following:

Table 1

The confusion table for defining classification statistics.

		Real labels
		Sensitive	Resistant
Predicted labels	Sensitive	TP	FP
	Resistant	FN	TN

It should be noted that the mentioned criteria are threshold dependent. AUPR and AUC are more reliable criteria that are independent of threshold values. AUPR is the area under the plot of Precision versus Recall at various thresholds. AUC is the area under the ROC Curve which plotted Recall against $$ FPR={\displaystyle \frac{FP}{FP+TN}}$$ at different threshold values.

The best threshold converting $$ \tilde{B}$$ estimated scores to binary labels were determined using Precision − Recall curve. The threshold value related to the elbow-point in this curve is considered as the best threshold because at this threshold, there is a good balance between Precision and Recall.

2.5 Variations of CDSML

The application of CDSML can be extended in order to be performed in various scenarios. The following subsections includes two variations of CDSML.

2.5.1 Performing on the response matrix with missing values

CDSML can handle the binary response matrix B without imputing missing values. To this aim, it is required to alter the Eq 12 such that it computes the loss function only for known pairs. Eq 21 is capable of handling the response matrix which contains missing values.

where Missing denotes the set of missing pairs. If we use Eq 21 instead of Eq 12, the method is able to be performed without applying the imputation step for missing IC50 values. The detailed formula for updating latent matrices in this version is explained in the S2 File. It should be noted that using these formula works also for response matrix without missing. For example, when the missing values are imputed, the set of Missing is empty; therefore, in that case the loss terms are calculated for all pairs.

2.5.2 Using double, single, or no similarity matrices

One can perform CDSML in three different scenarios based on similarity usage:

Double similarity: using both SC and SD similarity matrices

Single similarity: using either SC or SD similarity matrices

No similarity: using no similarity matrix.

To this aim, the loss function in Eq 21 must be modified. If only SC is ignored, the loss function will be changed to the Eq 22.

If only SD is ignored, the loss function will be changed to the Eq 23.

If both SC, SD are ignored, the loss function will be changed to the Eq 24.

It is notable that, in each scenario, the equations for updating latent matrices will be adjusted based on the related loss functions. The related equations for updating latent matrices in each scenario is presented in the S2 File.

It is notewothy that both variations explained in Sections 2.5.1 and 2.5.2 can be performed simultaneuosly. In other words, when the user wants to perform CDSML in each of similarity scenarios, the response matrix may contain missing values or not. Because the formula used in various similarity scenarios can ignore missing values if there are any.

3.1 Tuning hyper-parameters

We tuned CDSML hyper-parameters using grid search on different values of hyper-parameters on GDSC dataset using SC_E, SD_S. Then we used the same hyper-parameters for CCLE or when using other similarities. We considered α, β ∈ {0.125, 0.25, 0.5, 1, 1.5, 2, 2.5, ⋯, 8} and K = k′ * min(m, n), where k′ ∈ {0.1, 0.2, ⋯, 0.9}. The best values for hyper-parameters was determined based on AUC criterion because this criterion is independent from threshold value and assess the model more extensively. The best hyper-parameters for CDSML was α = 3.5, β = 4.5, and k′ = 0.7.

3.2 The performance of CDSML

CDSML predicts drug sensitivity according to the cell line similarity and drug similarities. We calculated three types of cell line similarities based on gene expression, mutation profile, and copy number variation. Each of these cell line similarities can be utilized as SC in CDSML similarity conservation terms. Furthermore, we computed three types of drug similarities based on chemical substructure, STiTCH network and target PPI, each of which can be considered as SD in CDSML similarity conservation terms. We analyzed the impact of using different types of similarities on the classification performance of CDSML. Table 2 represents the CDSML performance on GDSC dataset in three scenario, namely, double similarity, single similarity, and no similarity.

Table 2

The best value of each criterion is shown in bold.

The performance of CDSML on GDSC dataset using different scenarios including double similarity, single similarity, and no similarity.

	SC	SD	AUC	AUPR	Accuarcy	F1-score	Precision	Recall
No sim	-	-	0.8542	0.8846	0.7812	0.8291	0.7706	0.8974
Single similarity	SCE	-	0.9147	0.9372	0.8359	0.8653	0.841	0.891
	SCM	-	0.9148	0.9373	0.8356	0.8652	0.8401	0.892
	SCV	-	0.9071	0.9353	0.8288	0.8636	0.8346	0.8948
	-	SDS	0.9071	0.9353	0.8288	0.8636	0.8346	0.8951
	-	SDN	0.9072	0.9354	0.83	0.8637	0.8392	0.8899
	-	SDP	0.9071	0.9353	0.8292	0.8638	0.8357	0.8943
Double similarity	SCE	SDS	0.9158	0.9373	0.8388	0.8714	0.8426	0.9026
	SCE	SDP	0.9157	0.9373	0.8354	0.8714	0.8423	0.903
	SCE	SDN	0.9157	0.9398	0.8388	0.8715	0.8422	0.9031
	SCM	SDS	0.9148	0.9373	0.8357	0.8652	0.8402	0.8918
	SCM	SDP	0.9147	0.9373	0.8354	0.8652	0.8392	0.8929
	SCM	SDN	0.9147	0.9373	0.8350	0.8652	0.8388	0.8949
	SCV	SDS	0.9157	0.9397	0.8388	0.8714	0.8426	0.9026
	SCV	SDN	0.9157	0.9398	0.8387	0.8714	0.8424	0.9028
	SCV	SDP	0.9147	0.9372	0.8355	0.8653	0.8390	0.8934

We discuss about each of three scenarios in the following.

Discussing about “No Similarity” scenario:

The results in the first row of Table 2 indicate that CDSML gained high performance even using no similarities. CDSML performance using no similarity confirmed that the matrix factorization used in CDSML is efficient itself without using other extra information. The CDSML results using no similarity outperforms most of state-of-the-art methods mentioned in the paper.

Discussing about “Single Similarity” scenario:

Based on this table, one can conclude that using a single similarity matrix makes about 6% improvements of all criteria compared to the case of no similarities. Moreover, the performance of CDSML on various types of SC or various types of SD leads to almost equal performance, i.e. the calculated criteria for all executions on a single similarity matrix are the same. For example, using only the copy number variation similarity of cell lines leads to AUC of 0.9071 and AUPR of 0.9353, while using only the PPI similarity of drugs also leads to the same values of AUC and AUPR. Therefore, the impact of all similarity matrix on the CDSML performance is almost equal. Thus, CDSML yielded highly accurate and robust results.

Discussing about “Double Similarity” scenario:

Moreover, by comparing the results in single and double similarity scenarios, it can be inferred that the performance of CDSML using both SC and SD make subtle improvement in comparison to the single similarity scenario.

The performance of CDSML on CCLE dataset using different scenarios including double similarity, single similarity, and no similarity are provided in S2 File. Further evaluations on CDSML were conducted using cell line similarity based on gene expression and chemical substructure similarity of drugs.

3.3 Evaluation of predicted and imputed values for missing pairs

In order to show the rationality of imputed values in response matrix using Eq 1, one can compare the predicted labels by CDSML for missing values (denoted by $$ Pre{d}_L^{\left(CDSML\right)})$$ with imputed labels using Eq 1 (denoted by $$ \left(Imput{e}_L^{(Eq.1)}\right)$$. x L shows that the vector consists of binary labels. $$ Pre{d}_L^{\left(CDSML\right)}$$ is computed by executing CDSML on the binary response matrix with missing values (not applying imputation procedure), while considering all known pairs as training samples and all missing pairs as test samples. On other hand, $$ Pre{d}_L^{\left(CDSML\right)}$$ is computed by imputing missing values using Eq 1 and then converting the imputed values $$ Imput{e}_C^{(Eq.1)}$$ (index C denotes that the vectors contain the continuous values) to the imputed labels $$ Imput{e}_L^{(Eq.1)}$$ by comparing imputed IC50 values with max concentration thresholds of drugs. Since both vectors are binary, computing Accuracy, F1 − score, Precision, and Recall in addition to JI, cosine similarity, and cross entropy are meaningful and can give us an extensive comparison of these two vectors. Accuracy, F1 − score, Precision, and Recall are computed as defined in Section 2.4. The JI, cosine similarity, and binary cross entropy of two vectors X, Y are defined in Eqs 25–27.

The computed metrics by considering $$ Pre{d}_L^{\left(CDSML\right)}$$ as the ground truth are shown in Table 3. According to the significantly low value of cross-entropy as well as the high value of other metrics, one can conclude that $$ Imput{e}_L^{(Eq.1)}$$ and $$ Pre{d}_L^{\left(CDSML\right)}$$ are considerably similar. Consequently, the imputed values for missing pairs are reasonable.

Table 3

Comparison of $$ Imput{e}_L^{(Eq.1)}$$ and $$ Pre{d}_L^{\left(CDSML\right)}$$ on GDSC dataset.

Accuracy	F1-score	Precision	Recall	JI	Cosine similarity	Cross entropy
0.718	0.764	0.9297	0.6484	0.6181	0.7764	0.0065

To further justify the rationality of imputed values, we compared the imputed IC50 values using Eq 1 ($$ Imput{e}_C^{(Eq.1)}$$) with the predicted IC50 values by an state-of-the-art method. The method proposed by Zhang et al. predicts IC50 value uses a dual layer network which is similar to the idea used in Eq 1 [25], while having some differences. Moreover, Zhang et al. have shown that the method predicts reliable IC50 values for missing pairs by providing biological evidence for the missing IC50 values of three MEK inhibitor drugs. Thus, it is interesting to compare the IC50 values imputed using Eq 1 ($$ Imput{e}_C^{(Eq.1)}$$) with the predicted IC50 values by Zhang et al. method($$ Pre{d}_C^{\left(Zhang\right)}$$). To compare these two continuous vectors, regression criteria such as Root Mean Square Error (RMSE), Normalized Root Mean Square Error (NRMSE), and Mean Absolute Error (MAE) can be computed. RMSE, NRMSE, and MAE for two vectors X, Y are defined in Eqs 28–30.

Table 4

Comparison of $$ Imput{e}_C^{(Eq.1)}$$ and $$ Pre{d}_C^{\left(Zhang\right)}$$ on GDSC dataset.

RMSE	NRMSE(Impute)	NRMSE(Zhang)	MAE
0.04369	0.0065	0.0081	1.45

An interesting idea is to use Zhang et al. predicted IC50 values $$ \left(Pre{d}_C^{\left(Zhang\right)}\right)$$ for filling missing values in response matrix. To this aim, we converted continuous values of $$ \left(Pre{d}_C^{\left(Zhang\right)}\right)$$ into $$ \left(Pre{d}_L^{\left(Zhang\right)}\right)$$ with binary entries by comparing the max concentration with $$ Pre{d}_C^{\left(Zhang\right)}$$. In other words, we used Zhang et al. predicted labels $$ \left(Pre{d}_L^{\left(Zhang\right)}\right)$$ instead of labels computed by Eq 1; i.e. $$ Imput{e}_L^{(Eq.1)}$$ for filling the missing values in the response matrix. We then performed CDSML manifold learning on the obtained response matrix. Let us denote this version as CDML-Zhang and compared its results with CDSML using a stratified 5-fold cross-validation. Table 5 shows the comparison between CDSML and CDSML-Zhang. It can be seen that the evaluation criteria computed for both versions are so close to each other. However, CDSML leads to better results. These comparisons additionally validates using Eq 1 for imputing missing entries of response matrix.

Table 5

The assessments were done by averaging 30 repetitions of stratified five-fold cross-validation. The highest value of each criterion is shown in bold.

Comparison of CDSML and CDSML-Zhang performance on GDSC dataset.

Method	AUC	AUPR	Accuracy	F1-score	Precision	Recall
CDSML	0.9157	0.9398	0.8388	0.8715	0.8422	0.9031
CDSML-Zhang	0.912	0.937	0.836	0.870	0.84	0.896

To sum up all validation scenarios in this section and the related sections in S2 File, one can conclude that the results of all validations confirm that the imputed values using Eq 1 are reasonable and leads to the improvement in results.

3.4 Comparisons with classification methods

To compare the predictive performance of CDSML with other state-of-the-art classification methods, we used available implementations for HNMPRD [3], RefDNN [4], and DSPLMF [5]. In order to have a fair comparison, we evaluated these methods using 30 repetitions of stratified five-fold cross-validation on cell line-drug pairsfor GDSC and CCLE. These methods cover a variety of classification methodology and labeling thresholds. The methodology of HNMPRD, RefDNN, and DSPLMF are based on information flow, deep neural network, and matrix factorization, respectively. Moreover, the thresholds used in HNMPRD, RefDNN, and DSPLMF to convert IC50 values to sensitive/resistance labels are mixed Gaussian distribution, C_max, and drug-wise median, respectively.

A comparison between CDSML and other state-of-the-art methods for category classification is represented in Table 6. Considering GDSC dataset, HNMPRD obtained high Recall, but low values in other criteria. RefDNN showed satisfying performance according to all and DSPLMF obtained reasonable results according to Recall, but not high quality results with respect to other criteria. CDSML outperforms other classification methods by achieving much higher values of all criteria than other methods. It improves the results of RefDNN by almost 2% in AUC, 7% in AUPR, 2% in Accuracy, and 9% in F1 − score.

Table 6

The assessments were done by averaging 30 repetitions of stratified five-fold cross-validation. The highest value of each criterion is shown in bold.

Comparison of CDSML’s performance with state-of-the-art classification methods’ performance on GDSC and CCLE.

Dataset	Method	AUC	AUPR	Accuracy	F1-score	Precision	Recall
GDSC	CDSML	0.9157	0.9398	0.8388	0.8715	0.8422	0.9031
GDSC	HNMPRD	0.5728	0.6121	0.6088	0.7494	0.6032	0.9894
GDSC	RefDNN	0.9013	0.8753	0.8219	0.7828	0.8114	0.7583
GDSC	DSPLMF	0.7350	0.7218	0.6407	0.7096	0.5947	0.8797
CCLE	CDSML	0.9514	0.977	0.8989	0.9201	0.9485	0.8934
CCLE	HNMPRD	0.5817	0.7228	0.6527	0.7889	0.6528	0.9967
CCLE	RefDNN	0.6951	0.8200	0.6796	0.7825	0.7109	0.8713
CCLE	DSPLMF	0.8148	0.6944	0.7405	0.6667	0.5906	0.7668

Considering CCLE dataset, HNMPRD again obtained high Recall, but low values in other criteria. RefDNN and DSPLMF revealed acceptable performance. Consequently, the CDSML performance is significantly higher than other state-of-the-art classification methods and improves the best values of AUC, AUPR, Accuracy and F1 − score by more than 10%.

To fully compare the CDSML performance with classification methods, we compared its results with six off-the-shelf classification methods covering diverse methodologies: Gaussian Naiive Bayes (GNB), logistic regression (LR), random forest (RF), multi-layer perception (MLP), adaptive boosting (ADA), and K-nearest neighbor (KNN). The implementations of all these methods were conducted using the Scikit-learn python package [43]. It should be noted that the feature vector for each pair of cell line c_i and drug d_j was constructed by concatenating the ith row of SC and jth column of SD. Moreover, the C_max was used to label the cell line-drug pairs. The best hyper-parameter values were specified using grid search. The set of evaluated hyper-parameters and the best values of hyper-parameters for all methods are presented in Table 7. The best values of hyper-parameters were tuned using grid search and considering AUC criterion.

Table 7

The evaluated hyper-parameters and the best values of hyper-parameters for off-the-shelf classification methods.

Method	Evaluated hyper-parameters	Best hyper-parameters
GNB	Variance smoothing:{10−12, 10−9, 10−6, 10−3, 10−1}	Variance smoothing: 0.1
LR	Regularization scale:{ 0.001, 0.1,1, 0.01, 10, 100}	Regularization scale: 1
	Stop tolerance:{10−6, 10−4, 10−2}	Stop tolerance: 10−6
RF	Criterion: {gini, entropy}	Criterion: entropy
	Number of trees:{10,50,100,500,1000}	Number of trees: 100
SVM	Kernel: {linear, poly, RBF, sigmoid, precomputed}	Kernel: linear
	Regularization parameter: {0.01,0.1,1,10,100}	Regularization parameter: 0.1
MLP	Hidden layer sizes: {(50,50,50), (50,100,50), (100,)}	Hidden layer sizes: (50,50,50)
	Activation function: {tanh, ReLU}	Activation function: ReLU
	Solver: {SGD, Adam}	Solver: Adam
	Learning rate: {Constant, Adaptive}	Learning rate: Adaptive
	Regularization term: {0.0001, 0.05}	Regularization term: 0.05
ADA	Number of estimators: {10,50,100,500,1000}	Number of estimators: 50
	Learning rate: {1,1.25,1.5,1.75,2}	Learning rate: 1.25
KNN	K: {3,5,7,9,11,13,15,17,19,21,23,25}	K: 19
CDSML	α, β ∈ {0.125, 0.25, 0.5, 1, 1.5, 2, 2.5, ⋯, 8}	α = 3.5, β = 4.5
	K = k′ * min(m, n), where k′ ∈ {0.1, 0.2, ⋯, 0.9}	k′ = 0.7

Table 8 provides the comparison between CDSML and off-the-shelf classification methods. On both GDSC and CCLE, all methods showed good performance in classification of anti-cancer drug sensitivity; however, the least and highest values of criteria belongs to GNB and RF, respectively. On top of them, CDSML achieved the most accurate results and outperforms other methods according to all criteria, except Precision. Nevertheless, its Precision is not too far from the best Precision.

Table 8

The assessments were done by averaging 30 repetitions of stratified five-fold cross-validation. The highest value of each criterion is shown in bold.

Comparison of CDSML’s performance with off-the-shelf methods’ performance on GDSC and CCLE.

Dataset	Method	AUC	AUPR	Accuracy	F1-score	Precision	Recall
GDSC	CDSML	0.9157	0.9398	0.8388	0.8715	0.8422	0.9031
GDSC	GNB	0.8662	0.8974	0.7792	0.8033	0.8695	0.7465
GDSC	LR	0.9037	0.934	0.8243	0.8524	0.865	0.8401
GDSC	RF	0.9163	0.9435	0.8376	0.8653	0.8666	0.8641
GDSC	SVM	0.9038	0.9343	0.8243	0.8522	0.8663	0.8387
GDSC	MLP	0.9044	0.9351	0.8226	0.8533	0.8529	0.8562
GDSC	Ada	0.9039	0.9346	0.8268	0.8553	0.8629	0.848
GDSC	KNN	0.9091	0.9371	0.8341	0.8623	0.8646	0.86
CCLE	CDSML	0.9514	0.977	0.8989	0.9201	0.9485	0.8934
CCLE	GNB	0.9111	0.9523	0.842	0.8779	0.8914	0.8675
CCLE	LR	0.9435	0.9651	0.8811	0.9189	0.948	0.8916
CCLE	RF	0.9494	0.9632	0.8884	0.9193	0.9463	0.8939
CCLE	SVM	0.9444	0.9694	0.8867	0.9192	0.9478	0.8924
CCLE	MLP	0.9376	0.9671	0.8744	0.9146	0.9624	0.872
CCLE	Ada	0.9475	0.9609	0.8878	0.9118	0.9332	0.8924
CCLE	KNN	0.936	0.9653	0.8788	0.9116	0.907	0.8843

It is noteworthy that the computed criteria for machine learning models are significantly higher than the computed criteria for state-of-the-art classification and regression methods. The tuned hyper-parameters for machine learning models turn them into efficient and potent models that outperform other state-of-the-art methods. These findings are in agreement with the similar findings in RefDNN paper [4].

Figs 5 and 6 demonstrates the ROC curve and Precision − Recall curve for all of the classification methods on GDSC dataset, respectively. As it is shown the AUC and AUPR values for CDSML are 0.9221 and 0.943, respectively which are superior to the AUC and AUPR of other methods.

Fig 5

The ROC curve of all classification methods on GDSC dataset.

Fig 6

The Precision-Recall curve of all classification methods on GDSC dataset.

3.5 Comparison with regression methods

The prediction power of CDSML was further compared to the power of four state-of-the-art regression models: SRMF [6], CaDRReS [7], CDCN [8], and ADRML [9]. The implementation of these methods were available. We applied C_max threshold on the predicted IC50 values by these methods and convert them to the classification models. The performance of these models are provided in Table 9. On GDSC datset, CaDRReS achieved reasonable results. Moreover, SRMF and CDCN showed weak performance. Meanwhile, CDSML performs better than these methods with regard to all criteria. The interesting part of this evaluation is the comparison between CDSML and ADRML. According to Table 9 CDSML significantly outperforms ADRML. Since CDSML uses binary response matrix, it can be inferred that using binary response matrix as the initial input of manifold learning leads to more reliable classification of cell line-drug pairs into sensitive/resistant categories.

Table 9

The assessments were done by averaging 30 repetitions of stratified five-fold cross-validation. The highest value of each criterion is shown in bold.

Comparison of CDSML’s performance with state-of-the-art regression methods’ performance on GDSC and CCLE.

Datset	Method	AUC	AUPR	Accuracy	F1-score	Precision	Recall
GDSC	CDSML	0.9157	0.9398	0.8388	0.8715	0.8422	0.9031
GDSC	SRMF	0.4452	0.5493	0.712	0.7922	0.6908	0.9285
GDSC	CaDRReS	0.500	0.5922	0.6962	0.7738	0.6912	0.8788
GDSC	CDCN	0.4276	0.5460	0.7632	0.8164	0.7538	0.8903
GDSC	ADRML	0.4077	0.5096	0.7501	0.8177	0.7189	0.9481
CCLE	CDSML	0.9514	0.977	0.8989	0.9201	0.9485	0.8934
CCLE	SRMF	0.4539	0.6266	0.6966	0.8005	0.6998	0.9351
CCLE	CaDRReS	0.4389	0.6732	0.7133	0.8138	0.7048	0.9625
CCLE	CDCN	0.4262	0.6204	0.8608	0.8824	0.9796	0.8029
CCLE	ADRML	0.4257	0.6229	0.5473	0.6703	0.6373	0.7071

3.6 Tissue-specific conditions

Until now, we evaluated the methods on the whole dataset, which contains 19 tissue types. Nevertheless, most oncological treatments are designed based on tissue types [44], suggesting that considering tissue type may have a large impact on drug response predictions [3, 45]. Hence, we conducted tissue-specific assessments on three major tissue types.

The number of cell lines in each tissue type is shown in Fig 7. The most major tissue types are lung NSCLC, orogenital system, and leukemia with 67, 60, and 44 cell lines, respectively. We considered these tissue types and developed the tissue-specific models. The used cell lines in the train and test set for tissue-specific models belong to the same tissue type.

Fig 7

The number of cell lines in each tissue type.

The number shown on slices is the number of cell lines in the related tissue type. Three major tissue types are the offset slices.

Figs 8–10 illustrates the predictive performance of classification methods on three major tissues. The ranking of methods performance is similar to the methods learned on the whole dataset. All methods’ performance decline slightly due to the reduction in sample size, because the reduction in sample size limits the predicted power of models [46]. Nevertheless, CDSML outperforms other methods in three tissue-specific scenarios. The AUC of CDSML were 0.9096, 0.9186, and 0.9106 on leukemia, urogenital system, and NSCLC tissues, receptively.

Fig 8

The performance of methods in on NSCLC tissue.

Fig 9

The performance of methods in on orogenital system tissue.

Fig 10

The performance of methods in on leukemia tissue.

3.7 Case studies

As it was mentioned, some cell line-drug pairs have missing IC50 values, which were imputed in the pre-processing step. In order to conduct case studies, we considered the predictions for missing pairs and investigated their predicted novel sensitive pairs. There B matrix had 7790 missing values which accounts for roughly 40% of all samples. The predicted sensitivity scores for these pairs were obtained and sorted. The list of top 2000 most sensitive and top 2000 most resistant cell line-drug pairs are provided in S7 and S8 Tables in S1 File, respectively.

The top 15 most sensitive pairs that had missing associations in the original dataset were considered for further analysis. Table 10 represents the list of top 15 ranked samples assigned as sensitive. Reliable literature and the latest version of GDSC database were probed to provide evidence for the novel sensitive associations. As it is shown in Table 10, all top 15 novel associations were verified as sensitive pairs in the final version of GDSC. In addition, there are numerous insights about these associations in the literature.

Table 10

Top 15 novel sensitive pairs and pieces of evidence for these pairs.

Rank	Cell line name	Drug name	Cell line tissue	Sensitivity Score	Literature evidence	GDSC verification
1	CW-2	Ponatinib	large intestine	0.9253	[47]	verified
2	ZR-75-30	VX-702	breast	0.9248	[48]	verified
3	YAPC	Ponatinib	pancreas	0.9062	[49]	verified
4	NCI-H1092	VX-702	SCLC	0.9017	[50, 54]	verified
5	NCI-H1092	Temsirolimus	SCLC	0.8993	[51]	verified
6	CW-2	Vinorelbine	SCLC	0.8929	[52]	verified
7	NCI-H1563	VX-70	NSCLC2	0.8823	[53]	verified
8	COR-L88	VX-702	SCLC	0.8818	[50, 54]	verified
9	SHP-77	VX-702	SCLC	0.8783	[50, 54]	verified
10	LB373-MEL-D	Lenalidomide	melanoma	0.8758	[55, 56, 57]	verified
11	COR-L88	Temsirolimus	SCLC	0.8753	[51]	verified
12	CP66-MEL	VX-702	melanoma	0.875	[58]	verified
13	CW-2	Epothilone B	large intestine	0.8738	[59]	verified
14	NCI-H1092	Docetaxel	SCLC	0.8648	[60, 61]	verified
15	SCC-9	Ponatinib	head and neck	0.8619	[47, 62]	verified

These associations mainly report sensitive cell lines for Ponatini, VX-7002, Temsirolimus, Lenalidomide, Vinorelbine, Epothilone B, Docetaxel, among which, the insights about sensitive associations of three drugs are described in the following.

Ponatinib is a tyrosine kinase inhibitor which hinders the activity of four FGFR [63]. A recent study have shown that inhibits the cell growth in cell lines of various tissues types such as colon cancer [47]. Researchers states that the multi kinase inhibitors such as Ponatinib have showed efficient activity in targeting pancreatic cancer cells [49]. In addition, its effectiveness in the treatment of cell lines with head and neck cancer have been evaluated in several studies [47, 62].

VX-702—a P38 mitogen-activated protein kinase inhibitors—have been developed for the treatment of the inflammation diseases [64]. P38β MAPK is highly expressed in lung tissues and P38 MAPK pathways are highly activated in SCLC and breast cell lines, which leads to tumorgenesis and metastasis [48, 50]. The oral treatment with VX-702 seems to be effective in diminishing the fibrosis in SCLC and breast [50]. HuangFu et al. have shown that the administration of VX-702 along with INFβ treatments from melanoma, stabilization of the cell response, and improvement in the treatment efficacy [58].

Lenalidomide which has tumoricidal and immunomodulatory roles, has been frequently studied for the treatment of malignant melanoma and leads to high efficiency in combination with Docarbazine [55–57].

The supportive pieces of evidence in literature and GDSC database verified that CDSML could efficiently predict drug response label associations for the cell-line drug pairs.