Improved genome-scale multi-target virtual screening via a novel collaborative filtering approach to cold-start problem

Conventional one-drug-one-gene approach has been of limited success in modern drug discovery. Polypharmacology, which focuses on searching for multi-targeted drugs to perturb disease-causing networks instead of designing selective ligands to target individual proteins, has emerged as a new drug discovery paradigm. Although many methods for single-target virtual screening have been developed to improve the efficiency of drug discovery, few of these algorithms are designed for polypharmacology. Here, we present a novel theoretical framework and a corresponding algorithm for genome-scale multi-target virtual screening based on the one-class collaborative filtering technique. Our method overcomes the sparseness of the protein-chemical interaction data by means of interaction matrix weighting and dual regularization from both chemicals and proteins. While the statistical foundation behind our method is general enough to encompass genome-wide drug off-target prediction, the program is specifically tailored to find protein targets for new chemicals with little to no available interaction data. We extensively evaluate our method using a number of the most widely accepted gene-specific and cross-gene family benchmarks and demonstrate that our method outperforms other state-of-the-art algorithms for predicting the interaction of new chemicals with multiple proteins. Thus, the proposed algorithm may provide a powerful tool for multi-target drug design.

(b) Inability to adequately normalize the objective function in order to properly rank numerous solutions (i.e., ligands constructed de novo for the methods in the first group or ligands extracted from the compound libraries for the methods from the second group).
Recent years have seen the development of knowledge-based methods for protein-ligand interactions [15][16][17] . These algorithms rely on statistical and mathematical procedures to build upon the existing knowledge stored in the databases of known interactions 18 . In attempt to come up with more efficient and more accurate algorithms, biomedical researchers are starting to incorporate a variety of techniques from many different and seemingly unrelated fields. Recommender systems, which are used in the movie industry to predict users' preferences for movies, are finding their ways into computational molecular biology and biomedical research. In particular, techniques such as collaborative filtering 19 , compressed sensing 20 , and low-rank matrix completion 21 have been successfully applied to discover novel protein-protein interactions 22 and to reconstruct gene regulatory networks 23 . However, most of these methods have only sub-optimal performance in predicting preferences of new items. A computational method able to find targets for compounds with no available interaction data would help overcome the inaccuracy and complexity of de novo ligand design and virtual screening.
In this paper we present COSINE (COldStartINtEractions) -a statistical framework and a corresponding computational method for multi-target virtual screening via the "one-class collaborative filtering" technique. Our program exploits existing knowledge and databases of known interactions as well as the sequence similarities between proteins and structural similarities between drug molecules to suggest potential targets for new chemicals. Among unique aspects of our work are position specific weights, impute values, and a novel weighted-profile procedure for improving target prediction for novel chemicals. The accuracy of COSINE is validated in blind benchmarks that utilize well-known and publicly available resources. Our data shows that COSINE clearly outperforms numerous state-of-the-art methods for the same problem in several different tests and with respect to different accuracy measures. The algorithm is freely available at http://bioinfo.cs.uni.edu/ COSINE.html.

Methods
In a typical recommender system, user rating is expressed using different scores (e.g. 1-5 scale used by Netflix's users to rate movies). However, the nature of available data for protein-chemical interactions is different. Often times, only "positive" data consisting of known and validated interactions is available but there is no straightforward way of distinguishing "negative" scores (no interactions) from the missing data. The underlying binary score system (1 for interacting pairs and 0 otherwise) necessitates a deviation of the computational models used for protein-chemical interactions from the classical recommender models. COSINE belongs to the category of one-class collaborative filtering methods 24,25 since it does not treat all missing data as negative data. The protein-chemical interactions are predicted using the "low-rank matrix factorization" technique. More formally, given an incomplete matrix R of observed interactions, with m rows, representing targets, and n columns, representing chemicals, our algorithm decomposes R into a product of two lower dimension matrices of dimensions m × r and r × n, r ≪ min(m, n). The component matrices correspond to proteins' and chemicals' latent preferences. The assumption is that the set of proteins (respectively, chemicals) under consideration can be divided into a relatively small number of subsets with different proteins from the same subset exhibiting the same preferences to chemicals. Our algorithm takes account of the fact that related proteins, such as those with similar amino-acid sequences or similar three dimensional structures, exhibit similar preferences to chemicals and vice versa (structurally similar chemicals show similar preferences to proteins). Statistical framework. COSINE is a dual-regularized, one-class collaborative filtering method 25 that can employ either logistic or linear factorization. Our method can be thought of as a multi-directional extension of some recently described matrix factorization techniques for making recommendations 26,27 . Specifically, let m and n represent the number of proteins and chemicals, respectively, and let R = (r i,j ) be a m × n matrix of protein-chemical interactions In protein-chemical interaction studies, the binary matrix R is typically incomplete. While each nonzero entry r i,j = 1 signifies a known interaction, the meaning of each zero entry r i,j = 0 is ambiguous in that there can be either no interaction between the target t i and the compound c j , or, alternatively, that an interaction exists but it has never been verified experimentally. Thus, the goal is to predict the missing entries (i.e., to reclassify all unknown entries of the matrix R).
Building upon the general low-rank matrix factorization framework, COSINE approximates the probability of each chemical interacting with each target by mapping both chemicals and proteins to a common latent space of reduced dimensionality. The assumption here is that the number of factors influencing protein-chemical interactions is relatively small or, more formally, that the matrix of protein-chemical interactions is of low rank and, therefore, that it can be written as the product FG T of two matrices F and G of dimensions m × r and n × r, respectively, where r ≪ min (m, n) represents the number of latent factors. While our program is capable of performing either linear or logistic factorization, in the rest of this paper we will focus on logistic factorization, because it allows for an elegant statistical treatment.
Scientific RepoRts | 6:38860 | DOI: 10.1038/srep38860 Following Steck 28 , we first consider the loss function: i j T where f i and g j denote the i th and j th row (latent vector) of the matrices F and G, respectively, w i,j are the position specific weights on interaction scores, q i,j are the so-called "imputation values" 25 , λ F , λ G are tunable parameters and 2 denotes the Frobenius norm. The regularization terms F 2 2 and G 2 2 are included to prevent over-fitting. COSINE extends several other methods for the protein-chemical interaction prediction 16,27 , in at least two directions. Namely, the algorithm allows not only for the imputation of interaction values but also for different weighting of the interaction entries. In fact, to the best of our knowledge, COSINE is the only method for protein-chemical interaction prediction that employs position-specific weight and imputation values.
To provide insight into the motivation behind our method, consider, for instance, an ambiguous case where r i,j = 0 but some new experimental evidence suggests that the chemical c j might interact with protein t i . We can utilize this new knowledge by setting q i,j = 1 while lowering the corresponding weight w i,j to account for any uncertainty in the imputed value. A more thorough justification of the objective function (2) is given below. For a less general case, we refer the reader to Johnson 26 and Liu et al. 27 .

Position specific weights and impute values.
To derive the function (2) analytically, let e i,j be the event that the compound c j interacts with the target t i . Assume that the probability distribution of e i,j is logistic. In other words, assume that the probability p i,j assigned to e i,j is given by Recall also that w i,j reflects the confidence in the entry r i,j of the interaction matrix R. More precisely, higher weights are assigned to protein-chemical pairs (t i , c j ) which are known to interact (r i,j = 1), while lower values of w i,j are given to pairs for which r i,j = 0. To put it differently, a high number of positive training examples corresponds to each interacting pair while a lower number of negative training examples corresponds to each non-interacting (or unknown) pair. Hence, the likelihood of r i,j + q i,j given f i and g j can be written as , , or, at the matrix level, As in Steck 28 , the probability p(F, G|R + Q) can be derived through the Bayesian inference Finally, we derive the loss function (2) by taking the negative logarithm of (6), while assuming the Gaussian distribution of the entries of F and G 26 . Thus, in contrast to linear loss function 25 , namely the logistic loss function (used by default in our method) has an explicit probabilistic interpretation.

Dual regularization from proteins and chemicals.
To increase the accuracy of protein-chemical interaction prediction, we further extend the loss function (2) to account for the fact that similar chemicals are likely to interact with similar targets. Formally, let M = (m i,j ) be the matrix of pair-wise target similarity scores, where each entry m i,j represents the similarity between the proteins t i and t j , and let N = (n i,j ) be the matrix of pair-wise compound similarity scores. The affinity of similar chemicals to bind similar proteins is accounted for by minimizing the protein homophily Incorporating the regularization terms (8) and (9) above into (2), and introducing two additional tunable parameters, λ M and λ N , our loss function becomes i j T Figure 1 provides a toy example illustrating various components of the loss function. In practice, the entries m i,j of the matrix M typically represent the sequence similarity of the primary structures of proteins t i and t j , as measured, for example, by the normalized Smith-Waterman alignment score or by the PSI-BLAST e-value 29 . Alternatively, the values m i,j can be chosen to represent the three-dimensional similarity of the proteins' tertiary structures. Similarly, each n i,j represents the similarity score for the compounds c i and c j , as measured, for instance, by the Tanimoto score 30 or by the similarity of c i ′ s and c j ′ s pharmacological profiles 15 .
Note that the partial derivatives of (10) can be written as where  represents the Hadamard product. There are several ways to minimize the loss function (10) [25][26][27] . Similar to Liu et al. 27 , COSINE uses the AdaGrad -an iterative gradient descent method 31 . Weighted-profile approach for virtual screening. The most challenging task in protein-chemical interaction prediction is known as the "cold-start problem". The goal is to predict interactions of chemicals (or targets) for which no interaction data is available. COSINE implements a modified version of the "weighted profile" method 32,33 in which the latent preferences for a new protein (the rows of F) are computed as the sum of the latent preferences for that protein (calculated by the iterative minimization procedure, described above) and the latent preferences of J most similar proteins (those with available interaction data). More specifically, we set the i th row of the matrix F for the new target t i to where f j is the j th row of F (representing the latent preferences of the target t j ), v is the weight parameter and m i,j is the similarity score of the targets t i and t j . The normalization factor SM is set to The latent preferences for new chemicals (rows of G) are computed in the same way, using the compound similarity scores n i,j , namely.
Algorithmic details. COSINE minimizes the loss function (10) twice. The first time around, all imputation values q i,j are set to zero. The initial weights are set to 6 if r i,j = 1 and 1 otherwise, to reflect our increased confidence in experimentally verified interactions and lesser confidence in values r i,j = 0 (absent or unknown interaction). In the second iteration of the algorithm, the weight (which might be interpreted as the confidence in the value) of r i,j is increased by one if the computed probability of interaction p i,j is either too small or too large (more details are given in the Supplementary Table S1). The imputation values q i,j are adjusted in such a way that each entry of the new input matrix of interactions (namely r i,j + q i,j ) is set to 1 if the probability of the interaction p i,j computed in the first step is high (Supplementary Table S1). Otherwise, it is set to max(r i,j , p i,j ). We take max(r i,j , p i,j ) rather than p i,j since our underlying assumption is that the true interactions have been experimentally verified and hence the nonzero values of r i,j should be taken account of in the second step.

Results
To validate the algorithm, we compared it to a number of different methods for the same problem, namely KBMF2K 34  Benchmark #1. We first tested the accuracy of our algorithm in the classic Yam08 benchmark designed by Yamanishi et al. 32 . In this benchmarking experiment, which uses two different accuracy measures (AUPR and AUC), each dataset consists of four classes of targets: Enzymes, Ion Channels, GPCR's and Nuclear Receptors (Supplementary Table S2). In order to compare COSINE directly to the methods previously tested in this benchmark (KBMF2K, WNN and WNN-GIP) we performed a 5-fold cross-validation on the set of chemicals. More specifically, for each protein class, the set of chemicals was split into 5 subsets of approximately equal size and each subset was taken in turn as a test set. As described in van Laarhoven and Marchiori 33 , the training was performed on the remaining 4 subsets. The summary of the methods' accuracies, as measured by the area under the Precision-Recall curve (AUPR) and the area under the ROC curve (AUC), is given in Table 1. As seen in this table, while WNN method compares favorably to COSINE in the Enzyme class test, COSINE outperforms all of its competitors in all other target classes, most of the time, significantly. Moreover, the average AUPR and AUC scores achieved by COSINE exceed the average accuracies achieved by any other method tested in this benchmark.

Benchmark #2. Some methods for protein-chemical interaction prediction have been tested in Yam08
benchmark that uses 10-fold instead of 5-fold cross validation. To compare COSINE with those algorithms we modified the testing procedure and (similar to Liu et al. 27 ) ran 5 rounds of 10-fold cross-validation on the set of chemicals. Our findings are summarized in Table 2. As seen in these tables, COSINES achieves the best overall result, as measured by AUPR and AUC metrics. In contrast to 5-fold cross validation experiment, our method outperforms WNN-GIP in the Enzyme class benchmark with respect to both measures, but achieves a slightly lower AUPR than NRLMF (0.346 vs. 0.358).
It is interesting to note that COSINE's closest competitor in this test, namely NRLMF, also employs logistic factorization. However, unlike COSINE, the NRLMF method sets the weights globally, uses no imputation values and employs a different weighted profile scheme for cold start predictions. A different comparison of the two methods, using a different test sets and a different accuracy measure is presented in the subsection Benchmark #5 below.
Benchmark #3. Our next benchmarking data set, Yam 15 , has been constructed from the previous one 32 by extracting only the data corresponding to the compounds with available pharmacological profiles (Supplementary Table S3). Consequently, this benchmark mandates that all methods submitted use the similarity scores between pharmacological profiles computed by Yamanishi et al. 15  Strictly speaking, the only two algorithms that have been tested previously in the Yam10 benchmark using cross-validation on chemicals are the Yamanishi's 2008 algorithm, and its improved version, based on similarity of compounds' pharmacological profiles. Cobanoglu et al. have submitted their probabilistic matrix factorization method to a similar test 16  Benchmark #4. We have also compared COSINE to methods previously tested in the leave-one-out cross validation experiment that uses the Top 1 predictions as the accuracy criterion. Following the protocol described in the DASPfind paper 40 , for each target set and each drug under consideration, we removed all of the drug's known interactions and tried to retrieve them as Top1 predictions. As shown in Table 3, COSINE retrieves more interactions as Top1 predictions than any other method submitted to this benchmark. Although we have not trained the parameters of COSINE for this benchmark (we used the default ones found to work the best in the previous tests) it is reasonable to believe that the superior performance of COSINE over the other three methods is due to the fact that our algorithm has been explicitly developed to predict targets for new drugs (cold start). In contrast, the other three methods are tailored to not only "cold start" but also to "off-target" predictions.
Benchmark #5. Lastly, we compared the performance of COSINE to selected methods in the extensive ZINC benchmark. To generate the ZINC test sets, the ZINC data 43 was filtered by IC50 ≤ 10 μ M. This process yielded 31735 unique chemical-protein associations for 3,500 proteins and 12,384 chemicals. Cell-based tests and proteins appearing in protein complexes were excluded as well as proteins with unavailable primary sequences. Protein sequences were taken from UniProt 44 . Protein-protein similarity scores were calculated using BLAST. The ability of different algorithms to "rediscover" interactions was measured by "hiding" (setting to zero) the corresponding entries in the protein-chemical interaction matrix. To perform "cold start" analysis on ZINC data, we identified a set of chemicals having only one known target. The resulting set was further divided based on two criteria: 1) the number of chemicals the target proteins are associated with, and 2) the maximum chemical-chemical similarity score for the chemical in the dataset, with 0.1 increments. Each set was further subdivided into two subsets of approximately equal size, the test set (Supplementary Table S5) and the training set.
To provide for a conceptually different test, the ZINC benchmark uses the True Positive Rate (Recall or Recovery) at the top r% ∈ .
. . . . r ( {0 5, 1 0, 1 5, 2 0, 2 5}) of predictions for each chemical as the benchmarking measure. The Recall (Recovery) is defined as Recall = TP/CP, where TP and CP represent the total number of true and condition positives, respectively. Since there is a total of 3,500 targets, the r% of predictions include (35 · r) th or higher ranked target for each chemical. For instance, the True Positive Rate (TPR) of 0.7 at the 35 th cutoff rank (top 1%) means that 70% of the total tested positive pairs were ranked 35 th or better for the tested chemicals. Using TPR in place of AUC allows us to assess the performance of COSINE from a different angle. In particular, it is informative to compare COSINE head-to-head to NRLMF again, since, according to the results by Liu et al. 27 , the accuracies of NRLFM are significantly higher than KBMF2K, CMF, and WNN-GIP. Aside from NRLMF, we also analyzed the accuracy of one of the most popular and most widely used method for the cold start problem, based on the Parzen-Rosenblatt window (PRW) approach 38 . PRW is a highly accurate chemical structure-based target prediction method that uses neither the information obtained from proteins nor from the interactome. Finally, we submitted to the ZINC test a version of the COSINE algorithm, called REMAP, which uses linear (instead of logistic) factorization and global (instead of position specific) weights. Comparison with the latter algorithm (which has been used by our group for drug "off-target" prediction) is particularly useful since it illustrates the contribution of novel features of COSINE to its overall accuracy. Figures 2 and 3 demonstrate the performance of different algorithms in the ZINC test, as measured by Recall at the top r% of predictions for various values of r. The significant performance advantage of COSINE over REMAP illustrates the benefits of using local weights, logistic (as opposed to linear) factorization and a weighted profile approach for novel drugs and novel targets.
Additional analyses. We studied how the number of iterations in the matrix factorization step influences the accuracy of our algorithm. In our experiments, the convergence is attained after about 50 iterations for smaller data-sets (such as Nuclear Receptors test set) while a larger number of iterations (100-600) is needed to achieve comparable accuracy on larger data sets (such as Enzymes or ZINC test). As seen in Fig. 4, for very large data sets, such as ZINC, the added value has a low diminishing return after ~500 iterations. Thus, we opted for a reasonable speed-accuracy tradeoff of 600 iterations. Increasing the number of iterations further renders the algorithm computationally infeasible. A proper adjustment for the number of iterations results in the runtime comparable to other methods (Supplementary Table S6).
We also studied how COSINE performs in less than ideal settings, for instance, as a function of noise due to invalid or insufficient interaction data. We recorded the AUC values obtained on four target classes (NR, GPCR, Ion Channels and Enzymes) as a function of missing interaction data and as a function of incorrect interaction data. As shown in Fig. 5, our method is able to compensate a significant fraction of incorrect or missing data, due to the "low-rank matrix completion" technique built into the algorithm. This technique assumes that drugs' preferences to targets are determined by a relatively small number of interaction patterns. Explicitly imposing the rank constraint in the loss function (as done in COSINE and some other matrix factorization methods) results in eliminating erroneous interactions, those that cannot be explained by the small dimensionality of the space of latent preferences.

Discussion and Conclusion
Historically rational drug design has been characterized through identifying a single disease associated target and discovering exquisitely selective drugs against that target. Unfortunately, this one-drug-one-gene approach has been of limited success. This failure is manifest in the current issues facing the drug industry with near empty pipelines and costly post-market withdrawals. New methodologies are called for. Polypharmacology, which focuses on searching for multi-targeted drugs to perturb disease-causing networks instead of designing selective ligands to target individual proteins, has emerged as a new drug discovery paradigm 45 .
Computational methods that can assist polypharmacology are of key importance in drug development. In-silico protein-chemical interaction prediction has proven useful in drug-repurposing (drug-repositioning), an area of drug discovery that aims to find new therapeutic indications for known, FDA approved drugs 46,47 . Drug repurposing and other rational and structure-based drug design approaches are getting increased attention in the    pharmaceutical industry as the cost of bringing a new drug to the market is approaching $1 billion 48 . A significant portion of the drug development cost is attributed to the inability of many candidate drug compounds to pass stages II and III of clinical trials, which is due to their insufficient efficacy and/or increased toxicity. Hence, the drug discovery pipeline can be made more efficient by taking advantage of a systematic, rational approach. This strategy assumes an automated prediction and analysis of interactions on a large scale, carried out by comparing large subsets of the proteome against a wide array of existing and candidate drug compounds.
Selected statistical techniques, including recommender systems, known as "low-rank matrix completion" and "collaborative filtering", have been successfully used to predict protein-protein interactions 49 and to identify the gene clusters from the microarray data 50 . However, to date, the use of these systems in predicting protein-chemical interactions has been limited, due to their limitations in ability to accurately predict interactions of new compounds and new targets.
We introduce a computational method for predicting protein-chemical interactions based on matrix factorization. Our method builds upon "collaborative filtering" -a widely used statistical technique for making recommendations to utilize existing knowledge stored in the databases of known interactions. By incorporating  the weighting and imputation of the interaction data, as well as the dual regularization from both chemicals and targets, COSINE is able to exceed accuracy of other state-of-the-art methods for the same problem.
Our algorithm integrates chemoinformatics (chemical structural similarity), bioinformatics (protein sequence similarity) and a drug-target network (in form of matrix completion) approaches. Utilizing chemical structural similarity has proven useful (and has been widely applied) in the drug discovery for single-target virtual screening. Incorporating protein sequence similarity has shown promises in the prediction of drug off-targets [51][52][53][54] . Our drug-target network approach, formulated as a matrix completion problem, has been successfully applied to recommender system, which improve the performance of off-target prediction, especially when the chemoinformatics method fails.
The publically available chemogenomics data, on which all of existing virtual screening methods are inherently based, is incomplete and noisy. The missing interaction data is predicted in COSINE by completing the input interaction matrix, while biased and noisy data is filtered out by selecting the objective function that minimizes the rank of the output matrix of predicted interactions.
We recognize that, even though the ROC and PR curves may give a global estimation of the false positive rate for a prediction in the certain rank given by existing virtual screening algorithms, they may be not adequate for a risk-sensitive drug discovery application. In addition, in bio-and chemo-informatics applications, non-nested CV model is known to bias the parameters to the data set. Thus, other approaches to assessing reliability for specific new cases (including the label permutation and/or nested CV approach) will be extremely useful. We have developed several methods, e.g. ENTS 55 and case-based reasoning 56,57 for this purpose. In our on-going work, we plan to integrate these methods into the COSINE algorithm.