==== Front BioinformaticsBioinformaticsbioinformaticsBioinformatics1367-48031367-4811Oxford University Press 10.1093/bioinformatics/bty132bty132Original PapersSystems BiologyKernelized rank learning for personalized drug recommendation He Xiao 12http://orcid.org/0000-0002-5811-8875Folkman Lukas 12Borgwardt Karsten 12Wren Jonathan Associate Editor1 Machine Learning and Computational Biology Lab, Department of Biosystems Science and Engineering, ETH Zurich, Basel, Switzerland2 Swiss Institute of Bioinformatics, Basel, SwitzerlandThese authors contributed equally to this work. To whom correspondence should be addressed. Email: lukas.folkman@bsse.ethz.ch or xiao.he@bsse.ethz.ch15 8 2018 08 3 2018 08 3 2018 34 16 2808 2816 07 7 2017 4 1 2018 7 3 2018 © The Author(s) 2018. Published by Oxford University Press.2018This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.comAbstract Motivation Large-scale screenings of cancer cell lines with detailed molecular profiles against libraries of pharmacological compounds are currently being performed in order to gain a better understanding of the genetic component of drug response and to enhance our ability to recommend therapies given a patient's molecular profile. These comprehensive screens differ from the clinical setting in which (i) medical records only contain the response of a patient to very few drugs, (ii) drugs are recommended by doctors based on their expert judgment and (iii) selecting the most promising therapy is often more important than accurately predicting the sensitivity to all potential drugs. Current regression models for drug sensitivity prediction fail to account for these three properties. Results We present a machine learning approach, named Kernelized Rank Learning (KRL), that ranks drugs based on their predicted effect per cell line (patient), circumventing the difficult problem of precisely predicting the sensitivity to the given drug. Our approach outperforms several state-of-the-art predictors in drug recommendation, particularly if the training dataset is sparse, and generalizes to patient data. Our work phrases personalized drug recommendation as a new type of machine learning problem with translational potential to the clinic. Availability and implementation The Python implementation of KRL and scripts for running our experiments are available at https://github.com/BorgwardtLab/Kernelized-Rank-Learning. Supplementary information Supplementary data are available at Bioinformatics online. SNSF StartingSignificant Pattern Mining ==== Body 1 Introduction One of the key goals of precision medicine is the ability to suggest personalized therapies to patients based on their molecular profiles. This approach is particularly interesting for cancer treatment given the heterogeneous nature of this disease, with somatic alterations in cancer genes shown to be the determinants of a patient's response to therapy (Chapman et al., 2011). For the development of targeted treatments, large collections of patients with recorded clinical outcomes and molecularly profiled tumor samples are needed. However, the cost of such a medical records knowledge bank is currently prohibitive except for a few encouraging examples (Gerstung et al., 2017). Therefore, pre-clinical biological models, such as cultured human cancer cell lines, are a relatively inexpensive alternative approach for finding biomarkers. Recently, several large-scale drug sensitivity screens of genomically profiled cell lines ranging different cancer (sub)types have been established (Barretina et al., 2012; Iorio et al., 2016; Seashore-Ludlow et al., 2015). One of the current challenges lies in building accurate predictive models and translating these models into the clinic. A number of regression models (Costello et al., 2014; Jang et al., 2014) have been proposed to predict drug sensitivity measured by the half maximal inhibitory concentration (IC50) or the area under the fitted dose response curve (AUC). These methods typically build a predictor for each drug individually (Aben et al., 2016; Barretina et al., 2012; Iorio et al., 2016), use multi-task learning (Ammad-ud-din et al., 2016; Gönen and Margolin, 2014) or combine genomic profiles with information about the drug's chemical structure (Ammad-ud-din et al., 2014; Menden et al., 2013; Zhang et al., 2015). However, there are three potential problems with these approaches. Firstly, medical records knowledge banks are much sparser than cell line panels, meaning that only a few therapies are recorded for each patient compared to a comprehensive screen of tens to hundreds of drugs for each cell line. Secondly, these therapies are carefully selected by experienced doctors based on their expert judgment. Finally, the regression-based methods are not directly optimized for the prediction of the clinically relevant case in which a clinician needs to know a few most-suited drugs for a given patient. Here, we phrase the personalized drug recommendation as a ranking problem (Section 2) and propose a method, named Kernelized Rank Learning (KRL), that is directly optimized for the clinically relevant scenario of personalized drug recommendation illustrated above. To this end, we extend the approach of Weimer et al. (2007) and Teo et al. (2010) in Section 3. The key difference between KRL and the available regression approaches is that for a given cell line (patient), KRL directly predicts the top ranking drugs rather than the exact sensitivity values for all potential drugs. Our experimental evaluation using a dataset of cancer cell lines (Section 4) shows that KRL outperforms competing approaches, particularly in the realistic scenario when only a few training examples per cell line are available. Next, using publicly available patient data and the gene expression homogenization pipeline of Geeleher et al. (2014, 2017), we show that even when KRL is trained on cell line models, it can produce plausible recommendations for cancer patients. Finally, we summarize our findings and suggest potential directions for future work in Section 5. 2 Problem statement We are given X∈ℝn×p, a set of n samples (cell lines or patients) where each of them is represented with a molecular profile (e.g. gene expression) with p variables. Furthermore, we are given Y∈ℝn×m, a set of drug response measurements to m distinct drugs. Importantly, the majority of values in Y might be missing since the response of each individual sample is known only for a few of all potential drugs. Given X and Y as described above, and a molecular profile of a new sample, not present in X, our goal is to rank drugs based on their efficacy so that the most effective drug is ranked first. This goal is thus different from predicting the exact value of the drug response (as done by drug sensitivity prediction methods). Here we focus on the relative ordering (ranking) of a few, most effective, drugs. From a machine learning perspective, drug recommendation can be framed into a problem of learning a weight matrix W∈ℝp×m and then recommending drugs based on the predicted ranking vector f=xW. Here x∈ℝp is the molecular profile of a new sample, and f∈ℝm is a ranking vector which entails the predicted ranking scores for each of the m distinct drugs. Finally, the ranking of the drugs is obtained as a permutation π which sorts the drugs in the decreasing order of f (thus, π1=i such that max⁡(f)=fi, that is, π1 gives the identity of the first recommended drug, the one with the largest ranking score in f). The weight matrix W can be viewed as containing weights of all molecular features towards each drug. Since we are interested in relative ordering of the drugs, W needs to be learned jointly across all drugs rather then for each drug individually. Learning W can be formulated as the following minimization problem: (1) min⁡WL(XW,Y)+λR(W) where L(XW,Y) is a loss function, R(W) is a regularization term with λ controlling the effect of the regularization (regularization is necessary to avoid over-fitting to the training set). 2.1 Normalized discounted cumulative gain (NDCG) The key to our approach lies in employing a loss function based on the Discounted Cumulative Gain (DCG) ranking metric (Järvelin and Kekäläinen, 2002) as L(XW,Y) in Equation (1). Given the predicted ranking vector f∈ℝm and the true drug response vector y∈ℝm, DCG@k evaluates how well the order of the top k recommendations in f agrees with the order of the most effective drugs in y: (2) DCG@k(f,y)=∑i=1k2yπi−1 log⁡2(1+i) where π sorts f in decreasing order. DCG@k can be intuitively interpreted as follows: the numerator gives rewards proportional to the true responses of the k recommended drugs, while the denominator penalizes the incorrect ordering of these k recommendations. The gains are maximized when the most effective k drugs are the top k recommended drugs, ranked in the correct order. The Normalized DCG (NDCG) is more practical: (3) NDCG@k(f,y)=DCG@k(π,y)DCG@k(σ,y) where π and σ sort f and y in decreasing order, respectively. Hence, NDCG@k ranges from 0 to 1 and is maximized (perfect recommendation) when πi=σi,1≤i≤k. Finally, we convert NDCG@k, a gain function, into a loss function by taking 1−NDCG@k(f,y), which allows us to rewrite Equation (1) and state the drug recommendation problem formally: (4) min⁡W∑i=1n(1−NDCG@k(XiW,Yi))+λ||W||Frob2 where Xi (the i-th row of X) is the molecular profile of the i-th cell line and Yi (the i-th row of Y) is the true response vector of the i-th cell line to the m drugs. ||W||Frob2 is the squared Frobenius norm of W for the purpose of regularization. 3 Materials and methods We stated the drug recommendation as minimization of 1−NDCG@k. Unfortunately, it turns out that 1−NDCG@k is a highly non-convex and non-smooth function which makes its direct optimization difficult. Therefore, we minimize its convex upper bound derived by Weimer et al. (2007) instead and extend the Bundle Method for regularized Risk Minimization (BMRM) proposed by Teo et al. (2010) for the joint optimization of W across all drugs. Finally, Equation (4) employs a linear predictor (XiW), which does not take feature interactions into account. To boost prediction performance by introducing non-linearity, we apply the well-known kernel trick (Cortes and Vapnik, 1995). We thus refer to our approach as Kernelized Rank Learning (KRL). 3.1 Convex upper bound of the loss function Based on the work by Tsochantaridis et al. (2005), which describes how to find convex upper bounds of non-convex optimization problems, Weimer et al. (2007) showed that the following loss function is a convex upper bound of 1−NDCG@k(f,y): (5) ℓ(f,y)=max⁡π[1−NDCG@k(f,y)+〈c,fπ−f〉] where c∈ℝm is a decreasing sequence, e.g. ci=(i+1)−0.25 as suggested by the authors. Finding the permutation π is a linear assignment problem which can be solved with the Hungarian Marriage algorithm (Kuhn, 1955). 3.2 Kernelized rank learning (KRL) Non-linear kernel methods have shown good performance in a drug sensitivity collaborative competition (Costello et al., 2014). Therefore, we decided to kernelize the objective function defined in Equation (4). Let W=XTU, where U∈ℝn×m, the objective is rewritten as: (6) min⁡U∑inℓ(XiXiTU,Yi)+λTr(UTXXTU) where ℓ is the loss function defined in Equation (5), and Tr is the trace of a square matrix. By replacing XXT with a kernel matrix K on X, we get the KRL's kernelized objective function: (7) min⁡U∑inℓ(KiU,Yi)+λTr(UTKU) The time complexity of evaluating KRL's objective function comes down to the complexity of the linear assignment algorithm, which is O(m˜3), where m˜ is the number of non-missing drug response measurements for the given sample, which is typically only a few. Moreover, the KRL loss function can be decomposed and evaluated in parallel across all samples. Finally, to optimize Equation (7), we need the subgradient of the loss function L(KU,Y)=∑inℓ(KiU,Yi) with regards to U: (8) ∂UL(KU,Y)=K∂FL where F=KU and (∂FL)i=∂Fiℓ(Fi,Yi)=cπ¯−1−c where π¯ is the permutation that maximizes ℓ(f,y), π¯−1 is the inverse permutation of π¯, i.e. (fπ¯)π¯−1=f, and c is defined as in Equation (5). Supplementary Algorithm S1 details how to calculate the subgradient ∂UL. 3.3 KRL optimization The objective function in Equation (7) is convex but due to the maximization problem in ℓ expensive to evaluate. The BMRM algorithm proposed by Teo et al. (2010) has been shown to work well on such problems. However, the authors introduced only a single-task version of BMRM, which means it cannot be directly applied to solve Equation (7). Therefore, we extended BMRM for KRL to allow for the following differences when compared to the single-task BMRM: (i) U is a matrix rather than a vector and (ii) the objective function includes regularization on a positive definite kernel matrix, Tr(UTKU), rather than an L2 norm of a weight vector. Algorithm 1. BMRM for Kernelized Rank Learning  Data: U[0], K, Y, λ, ϵ>0, t=0  Result: U[t] 1 repeat 2   t=t+1; 3   A[t]=∂UL(KU[t−1],Y); 4   B[t]=L(KU[t−1],Y)−〈U[t−1],A[t]〉; 5   J(U)=λTr(UTKU)+max⁡10.05) to the second best method, KBMTL. Comparing KRL, which employs the RBF kernel, with the linear kernel implementation of our method (LKRL), KRL either outperformed (P<0.05, Precision@k for k∈{1,5,15,20}) or performed comparably to LKRL. Overall, KRL's performance was in the range of 23–36% for Precision@k and 47–58% for NDCG@k. Fig. 1. Comparison of KRL with related work (see Section 3.6 for details) in terms of Precision@k using the full training dataset for different values of the evaluation parameter k, which controls the number of predicted recommendations that are compared with the true drug ranking. The error bars show standard deviations from three cross-validation folds 4.2 Prediction using sparse training datasets To evaluate the compared methods in a more clinically relevant scenario, in which only a few therapies can be recorded for each patient, we subsampled the training dataset (training folds) at sampling rates of 50%, 20% and 10%. At the same time, we kept the test set (test folds) unchanged. In the case of 10% sampling, this strategy resulted in training on 22 drug responses per cell line on average. Figure 2 and Supplementary Figure S2 show Precision@5 and NDCG@5 as a function of the relative training dataset size, respectively. As expected, the performance of all compared methods decreased as the training dataset became sparser. Nonetheless, the Precision@5 improvements of KRL compared to the second method, KBMTL, were 1.1% (P=8.3×10−6), 2.1% (P=2.5×10−8) and 1.7% (P=5.2×10−7) for the sampling rates of 50%, 20% and 10%, respectively, compared to 0.6% (P=0.138) for the full training dataset. Thus, KRL appears to be more robust to sparse data than the regression approach. The same trend could also be observed for NDCG@5. The Precision@5 improvements of KRL compared to the other three methods using only 10% of the training dataset were 2.7% (P=2.5×10−7), 1.9% (P=2.4×10−7), 2.7% (P=4.9×10−8) and 5.5% (P=1.6×10−9) for LKRL, KRR, RF and EN, respectively. Fig. 2. Comparison of KRL with related work (see Section 3.6 for details) in terms of Precision@5 using the subsampled training datasets. The error bars show standard deviations from 10 randomly subsampled training datasets 4.3 Prediction using sparse training datasets biased towards effective therapies As discussed in the previous section, KRL appears to be more robust to missing data than related work. This is shown in Figure 2 where random subsampling of the training dataset was used to induce sparsity. Nonetheless, it can be argued that choosing a subset of drugs for each cell line randomly does not simulate medical records knowledge banks realistically. This is because clinicians do not prescribe therapies randomly but rather based on their expert judgment. Thus we employed another subsampling strategy: first, we selected a predefined proportion r of the most effective drugs per cell line and second, we randomly sampled q drugs per cell line from this selection. After this procedure, there might not be enough training data left for some of the drugs, therefore, after subsampling, we removed drugs with < 6 cell lines. We evaluated KRL and related work using training datasets generated with q∈{3,5,10} (yielding datasets with 143–199, 202–246, and 239–265 drugs in total, respectively) and r∈{100%,50%,20%,10%} (where r=100% is equivalent to sampling q drugs randomly). Here we analyze results with q=5 for brevity. Results for q=3 and q=10 showed a similar trend and are summarized in Supplementary Figures S3 and S4. Figure 3 and Supplementary Figure S3a show Precision@5 and NDCG@5, respectively, as a function of the fraction of the most effective drugs (r) used for sampling five drugs (q) per cell line. We found that as the training dataset was composed of more effective drugs, the improvement of KRL compared to the second best method, KBMTL, was gradually increasing from 0.7% to 6.2% (Precision@5) and from 0.3% to 10.4% (NDCG@5). Fig. 3 Comparison of KRL with related work (see Section 3.6 for details) in terms of Precision@5 using the subsampled training datasets, keeping five drugs (q) per cell line sampled from a predefined fraction (r) of the cell line's most effective drugs. The error bars show standard deviations from 10 randomly subsampled training datasets Both Precision@k and NDCG@k improvements of KRL compared to related work were statistically significant (P<0.05) for r∈{50%,20%,10%} across a wide range of k∈{3,5,10,15,20}. KRL was also able to outperform LKRL (P<0.05) across all evaluated values of r and k with the exception of r=100%, and r=10% for k = 1, in which case KRL and LKRL performed comparably. Figure 4 and Supplementary Figure S4a show Precision@k and NDCG@k as a function of the evaluation parameter k for the case of sampling five drugs from the 20% of the most effective drugs. Fig. 4. Comparison of KRL with related work (see Section 3.6 for details) in terms of Precision@k (for different values of the evaluation parameter k) using the subsampled training datasets, keeping five drugs (q) per cell line sampled from the 20% of the cell line's most effective drugs (r). The error bars show standard deviations from 10 randomly subsampled training datasets 4.4 Evaluation of the single top drug recommendations We found that KRL outperformed related work in terms of both Precision@k and NDCG@k for a variety of settings. Arguably, the most clinically relevant setting was presented in the previous section where a few (q) drugs were sampled from a given proportion (r) of the most effective drugs. We analyzed this set of results further by looking at each cell line's single top recommendation (i.e. the drug with the highest predicted ranking score) and quantifying how close this recommendation was to the most effective drug (i.e. the ground truth). To this end, we looked at the distribution of percentile ranks of the recommended drugs. We calculated a percentile rank of a recommended drug as the percentage of drugs to which the given cell line was less sensitive than it was to the recommended drug. For instance, a percentile rank of 90 means that the recommended drug was in the top 10% of the most effective drugs for the given cell line. Figure 5a–d compares histograms of percentile ranks for KRL and the second best method, KBMTL, when trained using five drugs (q) per cell line sampled from the 100%, 50%, 20% and 10% of the most effective drugs (r), respectively. The figure highlights that regardless of the sampling strategy, 49–56% of KRL's recommendations were within the top 10% of drugs to which the given cell line was most sensitive. While KRL and KBMTL performed comparably for sampling from the full dataset (Fig. 5a, P=0.50), the proportion of KRL's recommendations with percentile ranks >90 was higher by 4.5%, 18.6% and 24.4% for sampling from the 50%, 20% and 10% of the most effective drugs, respectively (Fig. 5b–d, P=7.3×10−5, 7.2×10−6 and 7.9×10−8, respectively). Supplementary Figure S5 shows a similar trend for sampling three and 10 drugs per cell line. Fig. 5. Histograms comparing the distributions of percentile ranks of drugs recommended by KRL and the second best method, KBMTL, using the subsampled training datasets, keeping five drugs (q) per cell line sampled from a predefined fraction (r) of the cell line's most effective drugs. The figure highlights that regardless of the sampling strategy, 49–56% of KRL's recommendations were within the top 10% of drugs to which the given cell line was most sensitive 4.5 Prediction with different data types So far, all abovementioned results were based on representing a cell line with its gene expression profile. It is of interest to evaluate KRL and related work also for the other three available molecular profiles: whole-exome sequencing, copy number variation and DNA methylation (see Supplementary Table S1 for a summary of features derived from these data types). Figure 6 and Supplementary Figure S6 compare Precision@5 and NDCG@5, respectively, for the six compared methods across all four available data types. These results are based on sampling five drugs (q) from the 100%, 50%, 20% and 10% of the most effective drugs (r) per cell line (Fig. 6a–d, respectively). In agreement with the findings from a drug sensitivity collaborative competition (Costello et al., 2014), gene expression was the most predictive data type regardless of the employed method. Focusing on the other three data types, KRL yielded a comparable performance to KBMTL (the second best method) for r=100% and 50% and achieved Precision@5 improvements of 3–9% (P<0.05) for r=20% and 10%. NDCG@5 improvements were in the range of 7–15%. Fig. 6. Comparison of KRL with related work (see Section 3.6 for details) across the four molecular data types: gene expression (GEX), whole-exome sequencing (WES), copy number variation (CNV) and DNA methylation (MET). The six compared methods were evaluated in terms of Precision@5 using the subsampled training datasets, keeping five drugs (q) per cell line sampled from a predefined fraction (r) of the cell line's most effective drugs. The error bars show standard deviations from 10 randomly subsampled training datasets 4.6 Clinical trials prediction with KRL Recently, Geeleher et al. (2014) showed that ridge regression models trained on in vitro cell line models can predict drug response in vivo. We were interested if KRL models can also generalize to patient data. We used the pipeline provided by Geeleher et al. (2014) to download and homogenize docetaxel breast cancer (Chang et al., 2003) and bortezomib multiple myleoma (Mulligan et al., 2007) clinical trials datasets. First, for docetaxel, there were 10 patients labeled as sensitive (S) and 14 as resistant (R). We used the predicted drug ranks to discriminate between the sensitive and resistant patients. We found that KRL had a tendency to rank docetaxel higher for the sensitive than the resistant patients (Wilcoxon rank sum test, P=0.04) with AUROC of 0.71 (Fig. 7a). Second, for bortezomib, there were 169 patients and their responses were classified as complete response (CR), partial response (PR), minimal response (MR), no change (NC), or progressive disease (PD). Comparably to the ridge regression approach (Geeleher et al., 2014), the medians of the predicted drug ranks were able to sort the five categorical responses in the correct order (Fig. 7b). The bortezomib trial defined also another classification which grouped 85 CR, PR and MR patients as responders and 84 NC and PD patients as non-responders. KRL was able to rank bortezomib higher for responders than for non-responders (Wilcoxon rank sum test, P=2.4×10−5) with AUROC of 0.68. Fig. 7. Predicted drug ranks for the docetaxel (a) and bortezomib (b) clinical trials. Docetaxel trial designated patients as sensitive (S) or resistant (R). Bortezomib trial classified the response into five categories: complete response (CR), partial response (PR), minimal response (MR), no change (NC) and progressive disease (PD). Moreover, the patients were grouped as responders (CR, PR and MR) and non-responders (NC and PD) 4.7 TCGA breast cancer cohort prediction with KRL We showed that drug ranks predicted with KRL trained on cell lines correlate with patient response to treatment. However, the clinical trials datasets do not allow for evaluating how well KRL ranks two (or more) alternative drugs for a given patient. Therefore, we employed the TCGA-BRCA cohort (The Cancer Genome Atlas Network, 2012) and analyzed different molecular breast cancer subtypes to gauge if KRL's drug recommendations for patients are plausible. Geeleher et al. (2017) showed how to homogenize gene expression profiles measured with RNA-seq (TCGA) and micro-arrays (GDSC), which in turn allowed them to train ridge regression models on the GDSC cell lines and use these models to ‘impute’ drug response of the TCGA patients. In their evaluation, they showed that HER2+ TCGA-BRCA patients were predicted to be more sensitive to lapatinib, a first-line therapy for HER2+ patients (Gomez et al., 2008), than patients without HER2 overexpression. Here we followed the idea of using molecular subtypes for evaluating KRL's recommendations and compared lapatinib to emerging targeted therapies for TNBCs (Kalimutho et al., 2015). First, we compared lapatinib and four PARP1/2 inhibitors (PARPi), veliparib, olaparib, talazoparib and rucaparib, some of which have shown promising therapeutic response in BRCA1/2-mutant (mBRCA) tumors (Tutt et al., 2010). Consistent with our hypothesis, we found that lapatinib was ranked higher than all four PARPi for 149 of the 163 HER2+ patients (91%). In contrast, lapatinib was recommended over (ranked higher than) PARPi only for one out of nine mBRCA TNBCs (11%). We re-trained KRL 10 times with randomly selected 100 drugs to confirm that the results were statistically different from a random recommendation rate of 50% (unpaired Student t-test, P=8.2×10−14 for HER2+ and 1.7×10−138 for mBRCA TNBCs). Table 1 details lapatinib recommendation rates for the individual PARPi. Upon randomly selecting 100 drugs, these results were also statistically different from random (P≪0.05) except for rucaparib (P=0.15), for which a recent clinical trial failed to report objective responses in breast cancer patients (Drew et al., 2016). Table 1. KRL’s recommendation rates for TCGA-BRCA molecular subtypes Recommendation HER2+ mBRCA TNBC JAK2+ TNBC n=163 n=9 n=56 lapatinib > PARPi 0.91 0.11 — lapatinib > veliparib 0.95 0.22 — lapatinib > olapariba 0.94 0.38 — lapatinib > talazoparib 0.92 0.22 — lapatinib > rucaparib 0.94 0.44 — lapatinib > ruxolitinib 0.86 — 0.29 Note: Recommendation rates shown here (lapatinib > alternative drug) are equal to the fractions of patients for which lapatinib was ranked higher than the alternative drug. a Olaparib was measured twice in the GDSC dataset but KRL performed comparably for both measurements, here we reported the mean recommendation rates. Second, we considered ruxolitinib, a JAK1/2 inhibitor approved for the treatment of myelofibrosis, currently being investigated in several clinical trials as a treatment for TNBC (Kalimutho et al., 2015). As expected, we found that lapatinib was ranked higher than ruxolitinib for 140 of the 163 HER2+ patients (86%), whereas it was recommended over ruxolitinib only for 16 of the 56 TNBC patients (29%) with normal-to-high JAK2 expression (JAK2+). We re-trained KRL 10 times with randomly selected 100 drugs to confirm that these results were statistically different from random (P=1.1×10−11 for HER2+ and 7.7×10−8 for JAK2+ TNBCs). 5 Conclusion and outlook This work phrased personalized drug recommendation as a ranking problem of choosing the most effective drugs. This follows from the observation that in a clinical setting, we are interested in providing a recommendation of a few, most effective, drugs rather than predicting the exact response (sensitivity) to all drugs available on the market. To this end, we proposed KRL, which directly optimizes a ranking loss function. In our empirical evaluation, we aimed to simulate a clinically relevant scenario, in which only a few therapies are prescribed for each patient, by subsampling the training cell lines dataset. We found that KRL outperforms state-of-the-art predictors in drug recommendation with the improvements being largest when the sampling was not random but biased towards effective therapies, thus simulating a medical records knowledge bank in which clinicians' prior knowledge is reflected. In order to provide reliable and robust recommendations in a clinical setting, KRL needs to be re-trained on patient data and treatment outcomes. Unfortunately, there is currently no publicly available dataset of sufficient size to evaluate KRL in this setting directly. Nevertheless, we were able to show that even when trained on cell lines, KRL recommendations can generalize to patient data and provide plausible recommendations. One potential direction for future work is extending KRL to learn treatment histories, interdependence and combinations of prescribed therapeutics. We envision that medical records knowledge banks which are currently being gathered by various research consortia will provide the opportunity to study these phenomena with a machine learning approach. Funding This work has been supported by the SNSF Starting Grant ‘Significant Pattern Mining’. Conflict of Interest: none declared. Supplementary Material Supplementary Data Click here for additional data file. ==== Refs References Aben N. et al (2016 ) TANDEM: a two-stage approach to maximize interpretability of drug response models based on multiple molecular data types . Bioinformatics , 32 , i413 –i420 .27587657 Ammad-ud-din M. et al (2014 ) Integrative and personalized QSAR analysis in cancer by kernelized bayesian matrix factorization . J. Chem. Inform. Model ., 54 , 2347 –2359 . Ammad-ud-din M. et al (2016 ) Drug response prediction by inferring pathway-response associations with kernelized bayesian matrix factorization . Bioinformatics , 32 , i455 –i463 .27587662 Barretina J. et al (2012 ) The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity . Nature , 483 , 603 –607 .22460905 Breiman L. (2001 ) Random forests . Mach. Learn ., 45 , 5 –32 . Broad Institute TCGA Genome Data Analysis Center (2016 ). Analysis-ready standardized TCGA data from Broad GDAC Firehose 2016_01_28 run. Chang J.C. et al (2003 ) Gene expression profiling for the prediction of therapeutic response to docetaxel in patients with breast cancer . Lancet , 362 , 362 –369 .12907009 Chapman P.B. et al (2011 ) Improved survival with vemurafenib in melanoma with BRAF V600E mutation . N. Engl. J. Med ., 364 , 2507 –2516 .21639808 Cortes C. , Vapnik V. (1995 ) Support-vector networks . Mach. Learn ., 20 , 273 –297 . Costello J.C. et al (2014 ) A community effort to assess and improve drug sensitivity prediction algorithms . Nat. Biotechnol ., 32 , 1202 –1212 .24880487 Drew Y. et al (2016 ) Phase 2 multicentre trial investigating intermittent and continuous dosing schedules of the poly(ADP-ribose) polymerase inhibitor rucaparib in germline BRCA mutation carriers with advanced ovarian and breast cancer . Br. J. Cancer , 114 , 723. 27002934 Garnett M.J. et al (2012 ) Systematic identification of genomic markers of drug sensitivity in cancer cells . Nature , 483 , 570 –575 .22460902 Geeleher P. et al (2014 ) Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines . Genome Biol ., 15 , R47. 24580837 Geeleher P. et al (2017 ) Discovering novel pharmacogenomic biomarkers by imputing drug response in cancer patients from large genomics studies . Genome Res ., 27 , 1743 –1751 .28847918 Gerstung M. et al (2017 ) Precision oncology for acute myeloid leukemia using a knowledge bank approach . Nat. Genet ., 49 , 332 –340 .28092685 Gomez H.L. et al (2008 ) Efficacy and safety of lapatinib as first-line therapy for ErbB2-amplified locally advanced or metastatic breast cancer . J. Clin. Oncol ., 26 , 2999 –3005 .18458039 Gönen M. , Margolin A.A. (2014 ) Drug susceptibility prediction against a panel of drugs using kernelized bayesian multitask learning . Bioinformatics , 30 , i556 –i563 .25161247 Iorio F. et al (2016 ) A landscape of pharmacogenomic interactions in cancer . Cell , 166 , 740 –754 .27397505 Jang I.S. et al (2014 ) Systematic assessment of analytical methods for drug sensitivity prediction from cancer cell line data. In Pacific Symposium on Biocomputing, pp. 63 –74 . Järvelin K. , Kekäläinen J. (2002 ) Cumulated gain-based evaluation of IR techniques . ACM Trans. Inform. Syst ., 20 , 422 –446 . Johnson W.E. et al (2007 ) Adjusting batch effects in microarray expression data using empirical Bayes methods . Biostatistics , 8 , 118 –127 .16632515 Kalimutho M. et al (2015 ) Targeted therapies for triple-negative breast cancer: combating a stubborn disease . Trends Pharmacol. Sci ., 36 , 822 –846 .26538316 Knijnenburg T.A. et al (2016 ) Logic models to predict continuous outputs based on binary inputs with an application to personalized cancer therapy . Sci. Rep ., 6 , 36812. 27876821 Kuhn H.W. (1955 ) The hungarian method for the assignment problem . Nav. Res. Logist ., 2 , 83 –97 . Leek,J.T. and Storey,J.D. (2007) Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genetics, 3, e161. Maxwell K.N. et al (2017 ) BRCA locus-specific loss of heterozygosity in germline BRCA1 and BRCA2 carriers . Nat. Commun ., 8 , 319. 28831036 Menden M.P. et al (2013 ) Machine learning prediction of cancer cell sensitivity to drugs based on genomic and chemical properties . PLoS One , 8 , e61318. 23646105 Mulligan G. et al (2007 ) Gene expression profiling and correlation with outcome in clinical trials of the proteasome inhibitor bortezomib . Blood , 109 , 3177 –3188 .17185464 Murphy K.P. (2012 ). Machine Learning: A Probabilistic Perspective . MIT Press , Cambridge, MA, USA . Pedregosa F. et al (2011 ) Scikit-learn: machine learning in Python . J. Mach. Learn. Res ., 12 , 2825 –2830 . Risso D. et al (2014 ) Normalization of RNA-seq data using factor analysis of control genes or samples . Nat. Biotechnol ., 32 , 896 –902 .25150836 Seashore-Ludlow B. et al (2015 ) Harnessing connectivity in a large-scale small-molecule sensitivity dataset . Cancer Discov ., 5 , 1210 –1223 .26482930 Teo C.H. et al (2010 ) Bundle methods for regularized risk minimization . J. Mach. Learn. Res ., 11 , 311 –365 . The Cancer Genome Atlas Network (2012 ) Comprehensive molecular portraits of human breast tumours . Nature , 490 , 61. 23000897 Tsochantaridis I. et al (2005 ) Large margin methods for structured and interdependent output variables . J. Mach. Learn. Res ., 6(Sep) , 1453 –1484 . Tutt A. et al (2010 ) Oral poly(ADP-ribose) polymerase inhibitor olaparib in patients with BRCA1 or BRCA2 mutations and advanced breast cancer: a proof-of-concept trial . Lancet , 376 , 235 –244 .20609467 Weimer M. et al (2007 ) COFI-RANK: maximum margin matrix factorization for collaborative ranking In Advances in Neural Information Processing Systems 20 , p. 1593 –1600 . Zhang N. et al (2015 ) Predicting anticancer drug responses using a dual-layer integrated cell line-drug network model . PLoS Comput. Biol ., 11 , e1004498. 26418249 Zou H. , Hastie T. (2005 ) Regularization and variable selection via the elastic net . J. R. Stat. Soc. B (Stat. Methodol.) , 67 , 301 –320 .