Drug cell line interaction prediction

TL;DR

This paper introduces tCNNS, a convolutional neural network model that predicts drug responses in cancer cell lines using drug structures in SMILES format, outperforming previous models.

Contribution

The paper presents a novel CNN-based model, tCNNS, that directly uses drug structures in SMILES format for phenotypic screening, achieving superior predictive performance.

Findings

Achieved R^2 of 0.84 and Pearson correlation of 0.92.

Outperformed previous models in drug response prediction.

Provided insights into phenotypic screening processes.

Abstract

Understanding the phenotypic drug response on cancer cell lines plays a vital rule in anti-cancer drug discovery and re-purposing. The Genomics of Drug Sensitivity in Cancer (GDSC) database provides open data for researchers in phenotypic screening to test their models and methods. Previously, most research in these areas starts from the fingerprints or features of drugs, instead of their structures. In this paper, we introduce a model for phenotypic screening, which is called twin Convolutional Neural Network for drugs in SMILES format (tCNNS). tCNNS is comprised of CNN input channels for drugs in SMILES format and cancer cell lines respectively. Our model achieves $0.84$ for the coefficient of determinant($R^2$) and $0.92$ for Pearson correlation($R_p$), which are significantly better than previous works\cite{ammad2014integrative,haider2015copula,menden2013machine}. Besides these statistical metrics, tCNNS also provides some insights into phenotypic screening.

Figures

Tables

Table 1. Table 1 : Tissue-Specific Test. The first column is the 13 tissue names and we range them in alphabetical order. The second column is the ground true IC 50 subscript IC 50 \text{IC}_{50} values for each tissue. The last two columns are the R 2 superscript 𝑅 2 R^{2} and R p subscript 𝑅 𝑝 R_{p} our model achieved by training using all the other tissue data. The number in the bracket is the result for the validation set.

Tissue name	Data Amount	$R^2$	$R_p$	RMSR
aero digestive	13806	0.703	0.843	0.0375
tract		(0.826)	(0.916)	(0.0280)
blood	31119	0.500	0.724	0.0449
		(0.833)	(0.917)	(0.0276)
bone	6826	0.659	0.813	0.0405
		(0.825)	(0.915)	(0.0283)
breast	9277	0.657	0.811	0.0383
		(0.829)	(0.919)	(0.0281)
digestive	17200	0.667	0.817	0.0384
system		(0.830)	(0.918)	(0.0282)
kidney	5199	0.669	0.819	0.0386
		(0.822)	(0.914)	(0.0286)
lung	34086	0.614	0.784	0.0371
		(0.827)	(0.919)	(0.0285)
nervous	15763	0.702	0.839	0.0364
system		(0.830)	(0.918)	(0.0280)
pancreas	5358	0.703	0.840	0.0370
		(0.820)	(0.913)	(0.0287)
skin	10488	0.676	0.824	0.0394
		(0.827)	(0.917)	(0.0281)
soft tissue	3165	0.712	0.853	0.0384
		(0.821)	(0.914)	(0.0284)
thyroid	2715	0.672	0.822	0.0410
		(0.833)	(0.918)	(0.0277)
urogenital	17112	0.715	0.849	0.0363
system		(0.825)	(0.914)	(0.0282)

Equations

$$y\mapsto \frac{1}{1+y^{-0.1}}.$$

$$S{S}_{\mathrm{tot}}=\sum _i(y_i-\bar{y}{)}^2$$

$$S{S}_{\mathrm{reg}}=\sum _i(f_i-\bar{y}{)}^2$$

$$S{S}_{\mathrm{res}}=\sum _i(y_i-f_i{)}^2=\sum _i{e}_i^2$$

$$R^2=1-\frac{S{S}_{\mathrm{res}}}{S{S}_{\mathrm{tot}}}.$$

$$R_p=\frac{\mathbb{E}[(Y-\bar{Y})(F-\bar{F})]}{{\sigma }_Y{\sigma }_F}.$$

$$\mathrm{RMSE}=\sqrt{\mathbb{E}[(Y-F{)}^2]}.$$

Full text

1 Introduction

Historically, drug discovery was phenotypic by nature. Small organic molecules exhibiting observable phenotypic activity (e.g. whole-cell activity) were detected, with a famous example being penicillin, which was serendipitously found. Phenotypic screening, an original drug screening paradigm, is now gaining new attention given the fact that in recent years the number of approved drugs discovered through phenotypic screens has exceeded those discovered through molecular target-based approaches. The latter, despite being the main drug discovery paradigm in the past 25 years, can potentially suffer from the failure in identifying and validating the therapeutic targets. In reality, most FDA approvals of first-in-class drugs actually originated from phenotypic screening before their precise mechanisms of actions or molecular targets ware elaborated. A popular example of this is aspirin (acetylsalicylic acid), for which it took nearly a century to elucidate the mechanism of its actions and molecular targets.

There are some public phenotypic screening datasets online to support the study the pharmacologic functions of drugs. Cancer Cell Line Encyclopedia (CCLE) and Genomics of Drug Sensitivity in Cancer (GDSC) are the most popular datasets in the field[4].

A pioneer work using machine-learning approaches to predict drug response on cancer cell lines was by Menden et al. [30]. The authors used a neural network to analyze the response of drugs to cancer cell lines on the GDSC dataset. Their main result was the achievement of $0.72$ for the coefficient of determination and $0.85$ for Pearson correlation. [3] and [16] are other two works on GDSC dataset. The first used kernelized Bayesian matrix factorization to conduct QSAR analysis on cancer cell lines and anti-cancer drugs, and the second used multivariate random forests. Both of their results are not as good as those in [30], which is chosen to be the benchmark for our work.

The first wave of applications of deep learning in pharmaceutical research has emerged in recent years. Its utility has gone beyond bioactivity predictions and has shown promise in addressing diverse problems in drug discovery. Examples cover bioactivity prediction[31], de novo molecular design[11], synthesis prediction[29] and biological image analysis[5, 26]. A typical example of applying deep learning in protein-ligand interaction is the investigation done by Ragoza et al[35].

Convolutional neural network (CNN) is a model in machine learning which can detect the patterns in data and support classification and regression[44]. CNN has achieved breaking-through results in many areas including pharmaceutical research [18, 40, 32] and has won the championship in ImageNet-2012[24].

Inspired by the achievements of CNN in these areas, we are interested to see if CNN, compared to conventional machine-learning techniques[30, 3, 16], could significantly improve the prediction accuracy of phenotypic drug response on cancer cell lines. CNN has not been applied to this problem. This is the first study in exploiting deep convolutional neural networks in predicting phenotypic anti-cancer drug response.

In this paper, we introduce a twin CNN networks model applied to drugs in SMILES format as input, and we call this model tCNNS. tCNNS comprises a CNN network for drugs and another CNN network for cancer cell lines, which will be explained in detail later. Comparing to previous work, we use the latest version of the GDSC dataset, which is bigger and more complete, and our model is more advanced. Most importantly, we achieve much better results than previous works. We share our model online, hoping to make a contribution to other researchers.

2 Relate Work

Recently, in drug discovery, researchers start to use the molecular structure of drugs directly as features [14, 13, 15, 2] instead of using extracted features from open source softwares[6, 38]. Due to their good ability to process high-dimensional structure on data, deep learning has been largely adopted in this area[25, 37, 9].

In [9], the author stated that both the qualitative classifiers and the quantitative structure-activity relationship (QSAR) models in this area depend on the molecular descriptor, which is the decisive step in the model development process. In the recent two years, there are several different deep neural network models that are trained directly from drugs structure and avoiding the decisive step. Those deep neural network models include unsupervised auto-encoder (AE), and supervised convolution neural network(CNN), recurrent neural network (RNN).

[14] and [20] converted simplified molecular input line entry specification (SMILES) of drug into vectors using unsupervised auto-encoder. Those vectors can be used as features or fingerprint of drugs. [2] and [41] further extended this method for drug discovery. They predicted the use of drugs by comparing the similarity between those vectors of drugs.

One of the direct ways to apply CNN on drug structures is to apply CNN on the image of drugs[12], instead of the formulas of drugs. In [12], the author adopts a computer vision method to screen the image of drugs. The advantage of starting from the image of drugs instead of the formulas of drugs is that it can avoid the massive work of handling the diversity of drugs. However, the disadvantages are that the accuracy is compromised because the information will be distorted when mapping drugs structures to images and the performance of this method relies on the quality of image processing.

Beside applying CNN on drug images, it is also possible to apply CNN on molecular 3D structures directly[39]. In this paper, the author predicted the binding energy of the small area around an atom, instead of the entire structure of drugs. We think that it will be interesting to compare the different representations of drugs, such as the 3D structured, the feature vectors learned from SMILES and the features extracted from other software like PaDEL [43]. They may have different influences on different problems. However, we did not find any kind of such comparisons.

Even though RNN is usually used to handle time sequence data[42] instead of spatial data, it is very impressive that [27] applied RNN on molecule SMILES to predict its solubility. They first converted the SMILES into indirect graphs, and then fed it into an RNN. In this work, they only considered the property of drugs alone, without considering the interactions among drugs and other biological factors, such as cell lines or proteins.

The benchmark work we used to compare with our model is [30], where the author used a neural network to analyze the $\text{IC}_{50}$ of drugs to cancer cells on the same dataset as ours. We think their network structure is not advanced enough, and the features they used are not informative enough. We design our tCNNS, a convolution neural network (CNN) based model, to predict the interaction between drugs and cell lines, and it performs better in the experiments.

3 Methods

In this section, we will introduce the chosen database GDSC with preprocessing steps and the adopted neural network structure in details to make our experiments easier to replicate.

3.1 Data Acquisition and Preprocessing

Genomics of Drugs Sensitivity in Cancer (GDSC) [8] is a public online database about the relationship among many types of cancer and various anti-cancer drugs. Cancer cell lines in GDSC are described by their genetic features, such as mutations state and copy number variances. For the drugs, GDSC provides their names and the compound id (CID). In chemistry, CID is a unique number assigned to each molecular and can be used as the reference number to extract more information about the drugs such as their molecular structures from other databases. GDSC uses $\text{IC}_{50}$ as the metric of drugs’ effect on cancers, which means how much of the drugs is needed to inhibit cancer by half. The less the value, the more effective the drug is. GDSC is an ongoing project and is being updated regularly. In this paper, we use version $6.0$ of GDSC, and the benchmark of our paper [30] used version $2.0$ of the GDSC.

We downloaded three files from GDSC:

1. (a)

Drug$\_$list.csv, which is a list of $265$ drugs. Each drug can be referred by its CID or name. 2. (b)

PANCANCER$\_$Genetic$\_$feature.csv, which is a list of $990$ cancer cell lines from $23$ different types of cancers. Each cell line is described by at most $735$ features. Any feature belongs to one of the two categories: mutation state or copy number alteration. 3. (c)

PANCANCER$\_$IC.csv, which contains the $\text{IC}_{50}$ information between $250$ drugs and $1074$ cell lines.

Note that the numbers of drugs in the file (a) and the file (c) are inconsistent, and the numbers of cell lines in the file (b) and the file (c) are also inconsistent. Some cell lines have less than $735$ features. Besides, GDSC does not provide the features for drugs, which have to be downloaded from other datasets. All of these indicate that a preprocessing is needed to clean the data.

The first step is to purify the drug list. There are $15$ repeating items in the file (a) which are removed. Some CIDs in the file (a) are inconsistent with the CIDs that in PubChem [23], which is a popular public chemical compounds database. To keep the consistency, we use the CIDs from PubChem. Some drugs cannot be found in PubChem by referring their names in the file (a) and they are removed. In the end, $223$ drugs with both names and CIDs are left.

The second step is to purify the cell lines list. For the $990$ cell lines in the file (b), $42$ of them has less than $735$ features. We remove these cell lines and there are $948$ cell lines left.

In the third step, we only keep the $\text{IC}_{50}$ values between the remaining drugs after the first step and the remaining cell lines after the second steps. All the other $\text{IC}_{50}$ values in the file (c) are removed. In summary, there are $223$ drugs and $948$ cell lines after the preprocessing. Amon the g $223\times 948=211,404$ interacting pairs, $81.4\%$ ($172,114$) of the $\text{IC}_{50}$ values are provided in file (c), and $18.6\%$ ($39,290$) of the $\text{IC}_{50}$ are missing.

The $\text{IC}_{50}$ data in file (c) are the logarithm of their real value. To make it easy for training and comparison, we follow the method in [30] to normalize the logarithmic $\text{IC}_{50}$ values into the $(0,1)$ interval. Given a logarithmic $\text{IC}_{50}$ value $x$, we first take exponential of it to get the real value $y=e^{x}$ and then use the following function to normalize $y$:

[TABLE]

Usually $y$ is very small ($<10^{-3}$), and the parameter $-0.1$ is chosen to distribute the result more uniformly on the interval $(0,1)$ [30].

3.2 Numerical Descriptor Extraction

Recently, there are some pioneering works that apply deep neural network (DNN) directly on the simplified molecular-input line-entry system (SMILES) of drugs. SMILES is a linear notation form to represent the structure of molecules, in which letters, digits and special characters are used to represent the chemical elements in a molecule. For example, ”C” is for carbon atom, ”=” is for covalent bond between two atoms, Carbon dioxide can be represent as O=C=O, and aspirin can be represented as O=C(C)OC1CCCCC1C(=O)O.

There are some challenges to apply convolution neural network (CNN) on drugs in SMILES format: first, SMILES can be constructed in various ways and there can be many possible SMILESs for each drug; second, the size of the samples for a CNN should be consistent, but the lengths the SMILES format of drugs are different from each other; third, and more importantly, the SMILES format is composed of different letters representing different chemical elements, such as atoms and bonds, and it does not make sense to apply convolution operation among different chemical elements. To solve these problems, preprocessing is needed to convert the SMILES into a uniform format, and restrict convolutional operation among the same chemical elements.

To keep unique SMILES format for the drugs, we use the canonical SMILES[34] for the drugs. Among $223$ drugs, the canonical SMILES of $184$ of them can be found from PubChemPy by the drugs names, which is a python interface for PubChem. The canonical SMILES of the remaining $39$ drugs are downloaded from the Library of the Integrated Network-based Cellular Signatures (LINCS) [21] instead.

The longest SMILES for the drugs contains $188$ symbols, and most SMILES lengths are between $20$ and $90$. To keep the size consistent and retain the complete information, all the SMILESs are left aligned with space padding in the right if they are shorter than $188$.

The neural network cannot directly take the drugs in SMILES format as input, and we need to convert the SMILES format (they are uniform length now after handling the second challenge) into a format that can be used in the neural network. There are $72$ different symbols in the SMILES format for the total $223$ drugs. The distribution of these symbols is quite imbalanced. For example, carbon atom [C] appears in all the $223$ drugs, but both the gold atom [Au] and chlorine atom [Cl] appear only once. Suppose we use rows to represent different symbols, and use columns to represent positions in the SMILES format, then each drug in SMILES format can be converted into a $72*188$ one-hot matrix which only contains [math] and $1$. The row number $72$ represents $72$ possible symbols, and column number $188$ represents the longest length of SMILES. In the one-hot matrix for a drug, a value $1$ at row $i$ and column $j$ means that the $i$th symbol appears at $j$th position in the SMILES format for the drug. In our model, the column of the one-hot matrix is treated as different channels in CNN, and the $1$D convolutional operation will be applied along each row of the one-hot matrix, which restricts convolutional operation within the same chemical element.

3.3 Deep Neural Network

The structure of our model is shown in Figure 1. The input data for our model is the one-hot representation of drugs (phenanthroline is used as an example for the drugs) and the feature vectors of the cell lines. The work-flow can be divided into two stages as follows.

During the first stage, we build a model with two convolutional neural networks (CNN) branches to distill features for drugs and cell lines separately. We use $1$D CNN for the cell-line branch since the input data is $1$D feature vectors for cell lines. We also use $1$D CNN for the drug branch and treat different symbols as different channels in CNN. The convolution is applied along the length of the SMILES format. The structures for the two branches are the same. For each branch, we use three similar layers: each layer with convolution width $7$, convolution stride $1$, max pooling width $3$, and pooling stride $3$. The only difference among the layers is that their number of channels is $40$, $80$ and $60$, respectively. The choices of those parameters for the CNN are inspired by the model in [22], in which the author chose a three-layers network model and used a prime number as the filter width. We find that either reducing the pooling size or adding the channel number has the potential to enhance our model but with the cost of losing stability. Losing stability means that the experiments results sometimes become un-repeatable and it will be explained in detail in Section 4.

In the second stage, after the two branches of CNN, there is a fully connected network (FCN), which aims to do the regression analysis between the output of the two branches and the $\text{IC}_{50}$ values. There are three hidden layers in the FCN, each with $1024$ neurons. We set the dropout probability to be $0.5$ for the FCN during the training phase [22].

We implemented our model on TensorFlow v$1.4.0$ [1], which is a popular DNN library with many successful applications [36, 1]. The data and code of this paper can be downloaded from https://github.com/Lowpassfilter/tCNNS-Project.

3.4 Performance Measures

We adopt the metrics of the coefficient of determination ($R^{2}$), Pearson correlation coefficient ($R_{p}$), and root mean square error ($\mathrm{RMSE}$) to measure the performance of our model, as in the benchmark paper [30].

$R^{2}$ measures variance proportion of the dependent variables that is predictable from the independent variables. Let $y_{i}$ be the label of a sample $x_{i}$, and our label prediction on $x_{i}$ is $f_{i}$. The error of our prediction, or residual, is defined as $e_{i}=y_{i}-f_{i}$. Let the mean of $y_{i}$ be $\bar{y}=\frac{1}{n}\sum_{i}y_{i}$, We have total sum of squares:

[TABLE]

regression sum of squares:

[TABLE]

residual sum of squares:

[TABLE]

$R^{2}$ is defined as:

[TABLE]

$R_{p}$ measures the linear correlation between two variables. We use $Y$ as the true label and $F$ as the corresponding prediction for any sample. Let the mean and standard deviation of $Y$ be $\bar{Y}$ and $\sigma_{Y}$ respectively, and those for the prediction $F$ be $\bar{F}$ and $\sigma_{F}$ respectively. $R_{p}$ is defined as:

[TABLE]

$\mathrm{RMSE}$ measures the difference between two variables $Y$ and $F$, and $\mathrm{RMSE}$ is defined as:

[TABLE]

4 Results and Discussion

In this section, we demonstrate the performance of our algorithm tCNNS with various data input settings, and it achieves high scores in the circumstances we set in the experiments.

4.1 Rediscovering Known Drug-Cell Line Responses

In the $223\times 948$ drug-cell line interaction pairs, GDSC provides the $\text{IC}_{50}$ for $172,114$ of them. In this part, we split those known pairs into $80\%$ as the training set, $10\%$ as the validation set, and $10\%$ as the testing set. In each epoch, parameters in tCNNS are updated using gradient descent on the training set. The validation set is used to check the performance of the tCNNS. If the $\mathrm{RMSE}$ on the validation set does not decrease in $10$ recent epochs, the training process stops and the predictions of our model on the testing set are compared with the given $\text{IC}_{50}$ values in GDSC.

The experiments are set in this way to stimulate the real situations in which we want to train a model on data with labels and apply the model on data without labels. The validation set is separated from the training set so that we can choose a suitable time to stop training independent and avoid the problem of over-fitting.

We check the regression accuracy of our model by comparing our predictions against the $\text{IC}_{50}$ provided in GDSC, and the results are displayed in Figure 2.

In the experiments, our model tCNNS outperforms previous work [30] in many ways: The efficiency of determination ($R^{2}$) is increased from $0.72$ to $0.84$, the Pearson correlation($R_{p}$) is increased from $0.85$ to $0.92$, and the root mean squared error ($\mathrm{RMSE}$) is reduced from $0.83$ to $0.027$.

Many hyper-parameters affect the performance of our model tCNNS, such as the number of layers and filter size. After examining, we find that smaller pooling size and more numbers of channels can enhance the performance furthermore but decrease the stability. Take the pooling size for example. When the it is reduced from $3$ to $2$, $R^{2}$ is further increased from $0.84$ to $0.92$ and $R_{p}$ is further increased from $0.93$ to $0.96$. The cost of this enhancement is that the network becomes unstable and diverge [19] during the training. To keep the experiments results repeatable, we report the results with parameters that ensure our experiments stable

In Figure 2, we notice that tCNNS is most accurate in the middle part, but less accurate in the two ends in the figure. Its outputs are not small enough when the experimental $\text{IC}_{50}$ values are very small and are not big enough when the data are huge. This means that tCNNS can be further optimized if we can enhance its performance in the two ends.

In the experiment, we do not recover the performance reported in [30] using the network structure in it. After replacing its network with tCNNS, We find that tCNNS do not converge using the features extracted from PaDEL. Finally, we replace the network in [30] with a deeper one, a network with three hidden layers and $1024$ neurons in each hidden layer. This modified benchmark gets the $R^{2}$ around $0.65$ and $R_{p}$ around $0.81$, which is shown on supplementary Figure 1. From this figure, we can see that the result is clearly horizontally stratified, which means that the neural network lacks representation power using PaDEL features.

Based on these results, we conclude that the connections among the $\text{IC}_{50}$ values and the SMILES format of drugs are stronger than that among the $\text{IC}_{50}$ values and features extracted using PaDEL.

4.2 Predicting Unknown Drug-Cell Line Responses

Previously, we split the known drug-cell line interaction pairs into training, validation, and testing set to check the accuracy of tCNNS. In this part, we predict the values for those missing pairs by training tCNNS on the known items. The known pairs are split into $90\%$ as the training set, and $10\%$ as the validation set. Again, if the $\mathrm{RMSE}$ on the validation set does not decrease in $10$ recent epochs, the training process stops and the trained network will be used to predict the values for the missing items. The results are shown in Figure 3.

Figure 3 is the box plot of the predicted $\text{IC}_{50}$ values for missing items grouped by drugs. For each drug, the box represents the distribution of the values with its related cell lines. The drugs are sorted by the median of the distribution: the $20$ drugs with highest median and $20$ drugs with the lowest median value are plotted. Since we do not know the real value for these missing pairs, the accuracy of our prediction is stated by survey and analysis as follows.

Bortezomib is the best drug in our prediction. In fact, the top $40$ pairs with the lowest $\text{IC}_{50}$ value are all from Bortezomib with some other cell lines. The outstanding performance of Bortesomib in missing pairs is consistent with that in the existing pairs. There is some supporting information in [7] that the author found that Bortezomib can sensitize cell lines to many other anti-cancer drugs.

Aica ribonucleotide and Phenformin have the poorest performance in our prediction. Based on our survey, the former one is initially invented to prevent blood flow, and the later one is initially used as an anti-diabetic drug. These two drugs have the potential to cure cancer because they can inhibit the growth of cell (Aica ribonucleotide) or inhibit the growth of Complex I (Phenformin), but their effects are limited since anti-cancer is only the side effect of these them, not their main function.

Based on our predictions, the $\text{IC}_{50}$ of drug Bortezomib with cell line NCI-H2342 is $1.19*10^{-4}\mu g$. This is the third smallest in our predictions and it indicates a good therapeutic effect. This prediction is supported by the findings in [10][17] that Bortezomib is able to control Phosphorylation that causes lung cancer and NCI-H2342 is a lung cell line. Similar evidence to supporting this prediction can also be found in Cell Signaling Technology’s 2011 published curation set. 111https://www.phosphosite.org/siteAction.action?id=3131

4.3 Retraining Without Extrapolated Activity Data

For each drug in GDSC, there are two important thresholds called minimum screening concentration (min_conc), which is the minimum $\text{IC}_{50}$ value verified by biological experiments, and maximum screening concentration (max_conc), which is the maximum $\text{IC}_{50}$ value verified by biological experiments. Any $\text{IC}_{50}$ beyond these two thresholds is extrapolated, and not verified by experiments. In general, $\text{IC}_{50}$ value within min_conc and max_conc is more accurate than those outside of the thresholds.

In the GDSC data that we use, only max_conc is provided, and there are $64,440$ $\text{IC}_{50}$ values below max_conc, which is about $37\%$ in the whole existing $172,114$ $\text{IC}_{50}$ values.

In this part, tCNNS is trained on the $\text{IC}_{50}$ values below the max_conc threshold. We keep $10\%$ data for validating and $10\%$ data for testing and reduce the percentage of data used for training from $80\%$ to $1\%$. The regression result is shown in supplementary Figure 2 and the comparison against the tCNNS which trained on whole existing data is showed in Figure 4.

From supplementary Figure 2, we see that tCNNS can achieve almost the same good result on max_conc data only, which is faster because less data is needed. There are some other properties of tCNNS that we can conclude from Figure 4. Firstly, it performs very well even with very limited training data. For example, even trained on $1\%$ of the existing $\text{IC}_{50}$ values, $R^{2}$ can be almost $0.5$ and $R_{p}$ be around $0.7$. Secondly, and more importantly, tCNNS performed better with less and more accurate data. We see that the dash lines (results on data below max$\_$conc) are always above the solid lines (result on all data), and the final performance on max$\_$conc data is almost as good as that on the total data, although the amount of data for the former is only $37\%$ of the later. To further compare the best performance on all data and max$\_$conc data only, we repeated the experiments for 20 times to see the distribution of the results, and show it in supplementary Figure 3.

There are three experimental results shown in Supplementary Figure LABEL:fig:max_conc, which are the experiment on all data, on the data below max$\_$conc and on a random subset of all data with the same size as those below max$\_$conc. Comparing the result on data below max$\_$conc with the result on the random data with the same size, we see that the performance of tCNNS is significantly better on data below max$\_$conc than on random data with the same size, and this proved that tCNNS is able to utilize the accurate data.

Comparing the result on data below max$\_$conc with that on all data, we see that the means of the tow $R^{2}$ are almost the same, and the $R_{p}$ on all data is only a little bit better than that on the data below mac$\_$conc. Meanwhile, the variations of $R^{2}$ and $R_{p}$ on data below max$\_$conc is a little bigger than those on all data. To conclude, the contribution from low quality extrapolated data is limited, and they can only reduce variation and improve $R_{p}$ a little bit.

Although the extrapolated data cannot enhance the accuracy of the model, they can help to improve the scalability. The tCNNS model trained on all data performs well when tested on data below max_conc, which is natural because the later is a subset of the former. However, the model trained on data below max_conc performs poorly when tested on all data. $R^{2}$ drops to $0.33$ and $R_{p}$ drops to $0.6$. Our explanation is that when tCNNS is trained on all data, it learns some general pattern among all data. On the contrary, when tCNNS is trained on data below max_conc, it only learns the specific patterns for this subset. So although the performance of tCNNS on all data and data below max_conc is similar, the paths they achieve the performance are different.

4.4 Blind Test For Drugs And Cell lines

In previous experiments, interaction pairs among drugs and cell lines are randomly selected to be in training set, validation set, or testing set, which means that a specific drug or a specific cell line can exist in training and testing at the same time, which is a problem of drugs re-positioning. In the real world, drugs discovery is a more challenging task in which we may not be able to train tCNNS on a specific drug or cell line before testing on them. Instead, we may need to predict the performance of the specific drug (or cell line) using tCNNS trained on other drugs(or other cell lines.) To stimulate this condition, we design the blind test for drugs and cell lines respectively.

In the blind test for drugs, drugs are constrained from existing in training and testing at the same time. We randomly select $10\%$ ($23/223$) drugs and keep their related $text{IC}_{50}$ values for testing. For the remaining $90\%$ drugs, we divide $90\%$ of their related $\text{IC}_{50}$ values for training and $10\%$ for validating.

In the blind test for cell lines, cell lines are prevented from existing in training set and testing set at the same time. Similar to the case of drugs, we randomly select $10\%$ ($94/948$) cell lines and keep their related $text{IC}_{50}$ values for testing. For the remaining $90\%$ ($904/948$) cells, we used $90\%$ of the related $\text{IC}_{50}$ for training and $10\%$ for testing.

We repeat the blind test for drugs on all data and on the data below max$\_$conc for $150$ times respectively to see the distribution of the results. The same number of experiments for the cell lines are also conducted. The results on all data are shown in Figure 5. The results on data below max_conc data are showed on supplementary Figure 4 respectively.

From Figure 5 and supplementary Figure 4 we see that the performance of tCNNS is comparably robust with the blind test for cell lines, but sensitive with the blind test for drugs. Without the knowledge of drugs in training, the performance drops significantly. Comparing the results in Figure 5 and supplementary Figure 4, we see that the extrapolated data makes no contribution in this setting.

To conclude, it is very important to have the specific drug or cell line in the training stage before we predict their performance. Experiments results support that even with only one or two related $\text{IC}_{50}$ value, the performance will be significantly improved. For example, NCI-H378 is a special cell line for lung cancer in GDSC that there are only two $\text{IC}_{50}$ values records for it. For other cell lines, all of them have at least $20$ $\text{IC}_{50}$ values. tCNNS can still make accurate predictions one of the $\text{IC}_{50}$ values for textitNCI-H378 if the other value is used during the training. Based our the results, the best drug for NCI-H378 are Bortezomib, which has been explained previously [10, 17].

Meanwhile, tcNNS predicts another potential drug Docetaxel for them. The predicted $\text{IC}_{50}$ value between Docetaxel and NCI-H378 is $0.03\mu g$ (third smallest for NCI-H378), and the predicted $\text{IC}_{50}$ between Docetaxel and NCI-H250 is $0.04\mu g$ (forth smallest for NCI-H250). It is reported in [33] that APR-246 is a potential useful drug on lung cancer because of its synergy with $TP53$ mutations in lung cells, and there is a ”additive effects” between APR-246 and Docetaxel.

Comparing the results of the blind test for drugs and blind test for cell lines, we guess the reason why the blind test for cell lines is slightly better is that there is more common information shared among different cell lines and less among drugs. For example, cell lines share similar genetic information, but drugs can be very diversified. To reduce the information sharing among cells lines, we designed another experiment in which cell lines from the same tissue cannot exist in training and testing at the same time. The result is shown in Table 1.

In GDSC, the $948$ cell lines belong to $13$ tissue types and $49$ sub-tissue types. We choose to use the $13$ tissue types instead of $49$ sub-tissue types because it can increase the distances and reduce the similarities among different tissues. Each time one tissue type is selected as testing data. For the rest tissues, we mix them together and use $90\%$ for training and $10\%$ data for validation. From Table 1, we can see that the performance decrease differently for different tissues. For example, blood has the lowest $R^{2}$ and $R_{p}$ in all tissues, which indicates that blood is the most different tissue from other tissues. From Table 1, we can also see that the result for validation is much better than testing, which not only because we use validating as the stop criterion for our model, but also because it shares more information with other tissues in the training stage.

The above results in the blind test give us some hint that with limited budgets, we should carefully arrange the in vivo experiment to cover a wider range of drugs and cells from different tissues to get better in silicon predicting power.

4.5 Cell Lines Features Impacts

In GDSC, the $735$ features for cell lines after preprocessing belongs to $310$ gene mutation states, and $425$ copy number variations. Since different laboratories may use different methods to extract the features for cell lines, in reality, it is not easy to have the complete $735$ features for all cells. More possibly, researchers may have smaller and different features group for cell lines. It will be attractive if our tCNNS can have good performance with fewer features for cell lines. In this part, we test tCNNS performance with a different number of features for cell lines. All the results in this part are shown in Figure 6.

4.6 Biological Meaning v.s Statistical Meaning

tCNNS takes the one-hot representation of the SMILES format as the features for drugs. Initially, in the one-hot representation of the SMILES format, each row represents a symbol, and each column represents a position in the SMILES format, which is left aligned. For human researchers, the SMILES format is a well-defined concept with biological meaning. However, we think tCNNS may lack the ability to comprehend the biological meaning of the SMILES format and it instead relies on the statistical pattern inside the data. To verify this hypothesis, we modify the one-hot representation of the SMILES format in three ways. Firstly, we randomly shuffle the order of the symbols, which equals shuffling the rows in the one-hot representation. Secondly, we cut the SMILES format into two pieces and switch their positions, which equals shifting the columns in the one-hot representation. Thirdly, we shuffle the positions in the SMILES format, which equals shuffling the columns in the one-hot representation. In the last two ways of the modification, the biological meaning of SMILES is corrupted. We repeat the experiments in the three settings for $10$ times respectively and compare the results to the result using the SMILES format without any modification as the benchmark. The comparison is shown in supplementary Figure 5.

Since the computation of different channels (the rows in the one-hot representation of SMILES format) are independent of each other in the convolution neural network, we think the result should keep the same in the first modification with the benchmark. In the second modification, shifting the SMILES will only change the three column in the one-hot representation: the two ends columns and the column where we start the shifting. Meanwhile, other columns keep the same. Thus, we also think that the result in the second modification should keep the same with the benchmark. The experimental results do support our predictions in the first two modifications. What surprised us is that the performance also keeps similar in the third modification. We used to think that shuffling the positions in the SMILES format will deteriorate the result because it corrupted the structure of the SMILES format. This stability means that tCNNS actually does not capture the biological meaning of the SMILES format for drugs, and it relies on the statistical patterns inside the SMILES format, cell line features, and the $\text{IC}_{50}$ values.

4.7 Eliminating Outlier

In the last column in supplementary Figure 5, we compared the results of tCNNS with the results of benchmark [30]. As GDSC has changed in recent years, we can not use the same data as [30]. In the experiment, we follow the method introduced in [30] and apply it to current data. We use PaDEL(version 2.1.1) to extract $778$ features for each drug. For cell lines, we use the same $735$ features, instead of the 157 features in the old version of GDSC used in [30].

To compare the features extracted using PaDEL and features extracted from the SMILES format using CNN, we visualize the distribution of the drugs in different feature spaces to check the difference. When drawing the distribution of drugs using CNN, we use the output of the last layer of CNN tranche for drugs. The reason is that in deep neural network, the fully connected layer is responsible for regression analysis, and CNN is for extracting the high-level features from the drugs features. The input data for the fully connected network is the output of CNN tranche. So drawing the distribution of drugs using the output of the last layer in CNN is more reasonable.

We also compared the distribution of cell lines in features space of GDSC, and in the output space of the last layer in CNN respectively. The visualization tool we used is t-SNE[28], which is widely used to visualize high dimensional data in deep learning.

The visualization results are shown in Figure 7. We can see that CNN can distribute the drugs and cell lines more uniformly than features using PaDEL and features of GDSC. For drugs, we see that there are seven outliers in the PaDEL features space. These seven outliers are the only seven drugs that are composed of multiple parts. We show the structure of those drugs in supplementary Figure 6.

5 Conclusion

In this paper, we came up with a model tCNNS for phenotypic screening between cancer cell lines and anti-cancer drugs. tCNN is tested on a new version of GDSC with more data compared to previous works. It achieved a much better coefficient of determinant and Pearson correlation than previous works and made predictions for missing values in GDSC with trustful evidence. tCNNS can also converge with a very small set of training data and fewer features for cancer cell lines, which is economically efficient. tCNNS took SMILES as the input data for drugs, and this eliminated the outlier problem in previous works in which drugs fingerprint are used as the features. tCNNS performs better in drugs repurposing than drug discovery, this may be a future research topic.

6 Competing interests

The authors declares that they have no competing interests.

7 Declaration

The authors declares that they have provided the code and data public accessible.