Mateo et al. Genome Medicine (2020) 12:78 https://doi.org/10.1186/s13073-020-00774-x METHOD Open Access Personalized cancer therapy prioritization based on driver alteration co-occurrence patterns Lidia Mateo1, Miquel Duran-Frigola1, Albert Gris-Oliver2, Marta Palafox2, Maurizio Scaltriti3,4, Pedram Razavi3,5, Sarat Chandarlapaty3,5, Joaquin Arribas6,7,8,9, Meritxell Bellet10,11, Violeta Serra2,9 and Patrick Aloy1,8* Abstract Identification of actionable genomic vulnerabilities is key to precision oncology. Utilizing a large-scale drug screening in patient-derived xenografts, we uncover driver gene alteration connections, derive driver co-occurrence (DCO) networks, and relate these to drug sensitivity. Our collection of 53 drug-response predictors attains an average balanced accuracy of 58% in a cross-validation setting, rising to 66% for a subset of high-confidence predictions. We experimentally validated 12 out of 14 predictions in mice and adapted our strategy to obtain drug- response models from patients’ progression-free survival data. Our strategy reveals links between oncogenic alterations, increasing the clinical impact of genomic profiling. Keywords: Drug-response biomarkers, Driver co-occurrence networks, Precision oncology Background genes published to date has analyzed over 9000 tumor In light of the complexity and molecular heterogeneity samples, across 33 tissues of origin, and has systematic- of tumors, clinical and histopathological evaluation of ally identified driver mutations in 258 genes [6]. cancer patients is nowadays complemented with gen- Approximately half (142) of those driver genes were as- omic information. Genome-guided therapy has been sociated with a single tumor type, whereas 87 genes shown to improve patient outcome [1, 2] and clinical seem to provide a growth advantage in several tumor trial success rate [3], and despite some controversy [4], types. The number of drivers detected per tumor type prospective molecular profiling of personal cancer ge- varies widely, ranging from 2 in kidney chromophobe nomes has enabled the identification of an increasing cancer to 55 in uterine cancer. Despite the large number number of actionable vulnerabilities [5]. of drivers identified per tumor type, every patient has a Cancer genome sequencing initiatives have found that unique combination of mutations and copy number vari- any given tumor contains from tens to thousands of mu- ants: 90% of patients show at least one putative driver tations. However, only a few of them confer a growth alteration, but each sample only contains a median of advantage to cancer cells, driving thus the tumorigenic three putatively altered drivers [7]. process. The most comprehensive study of “driver” On top of identifying key alterations in tumor develop- ment, it is fundamental to pinpoint those that can shed * Correspondence: patrick.aloy@irbbarcelona.org 1 light on the most appropriate therapy to treat each Joint IRB-BSC-CRG Program in Computational Biology, Institute for Research tumor (i.e., biomarkers). Often, patients with similar in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology, Barcelona, Catalonia, Spain clinicopathological characteristics might be molecularly 8 Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, different [6]; this inter-patient heterogeneity is one of Catalonia, Spain the reasons why only a subset of them will actually Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Mateo et al. Genome Medicine (2020) 12:78 Page 2 of 23 respond to a given targeted treatment. Computational MAP2K1 and upstream activating mutations in BRAF or studies suggest that up to 90% of patients may benefit NRAS) [16]. Additionally, at pathway level, the concomi- from molecularly guided therapy when biomarkers of tant activation of PI3K signaling pathway with FGF sig- uncertain clinical significance, as well as off-label and naling (FGFR2 and FGFR3) or with NRF2-mediated experimental drugs, are used to guide treatment selec- oxidative response has also been identified in several tion [7, 8]. Although randomized controlled trials are tumor types [16]. In this context, a single gene-based still considered the gold standard in the clinics, they stratification of patients into subtypes and treatment cannot address all possible patient clinicopathologic and arms might be over-simplistic, and novel frameworks molecular subtypes [9]. Precision medicine has that exploit co-mutational patterns might prove more prompted the reconsideration of clinical drug develop- effective. ment pipelines, with the implementation of more sophis- As in the identification of driver mutations, the discov- ticated clinical trial designs, such as umbrella, basket, ery of drug-response biomarkers requires large numbers and platform trials to account for inter-patient hetero- of patient molecular profiles matched to treatment out- geneity [10]. In particular, the implementation of adap- comes. Unfortunately, treatment history information of tive enrichment strategies allows for continual learning large-scale genomics endeavors has not been systematic- and modification of the eligibility criteria as data accu- ally collected (e.g., TCGA [23]) or is not yet publicly mulate, with the objective of recruiting those patients available (e.g., GENIE Consortium [24]). Even though that are most likely to benefit from treatment [9–12]. better data sharing policies are needed, many concerns However, despite the implementation of these novel ex- are raised regarding privacy, property, and the prelimin- perimental designs, currently, only alterations in 28 ary nature of confidential biomedical data. Safer alterna- genes have accumulated enough clinical evidence to be tive ways of sharing biomedical data are already on the approved as biomarkers by the FDA [13]. Indeed, a re- table [25], but until the access to systematically anno- cent comprehensive analysis of 6729 pan-cancer tumors tated clinical records becomes a reality, the research could only identify actionable mutations with thera- community largely relies on drug-response data gathered peutic options available in clinical practice (FDA-ap- from pre-clinical models. proved or international guidelines) or reported in late Cancer cell lines are the most widely used in vitro phase (III–IV) clinical trials in 5.2% and 3.5% of the model system and have been fundamental tools to set samples, respectively [14]. These figures coincide with the grounds of our understanding of cancer biology and clinical trial enrolment rates [1], where only 89 out of to assess the efficacy of a broad spectrum of cancer 1640 of patients could enter genotype-matched treat- drugs [26]. Unfortunately, cancer cell lines have been ment trials, the vast majority of which involved muta- cultured as monolayers on plastic surfaces, and in tions in four genes, namely PIK3CA, KRAS, BRAF, and growth-promoting conditions, for decades. As a conse- EGFR. This highlights an acute need to expand the quence, most of them have suffered a substantial tran- current repertoire of response biomarkers to cover more scriptional drift, and they likely represent a cell drugs and patients. subpopulation from the original primary tumor [27]. The eligibility criteria of most genomically matched Those facts have fueled the debate regarding how well basket clinical trials are based on the single-gene bio- cancer cell lines resemble the tumors from which they markers. However, most tumors do not present a single were established and to which extent they are clinically actionable mutation but have co-occurring driver alter- relevant [15, 27]. A more realistic model to bridge the ations that might simultaneously alter key players of sig- bench-to-bedside gap is the patient-derived mouse xeno- naling pathways connected by cross-talk and feedback graft (PDX) [28]. To some extent, PDXs preserve inter- mechanisms [15, 16]. There are many documented cases and intra-tumoral heterogeneity, and mimic the clinical of functionally relevant co-occurring oncogenic muta- course of the disease and response to targeted therapy, tions, such as the concomitant inactivation of TP53 and at least in certain tumor types [29–31]. Indeed, a recent RB1 [17], co-deletion of CDKN2A and CDKN2B [18], review reported a 91% (153 out of 167) correspondence co-amplification of MDM2 and CDK4 [19, 20], 1p/19q between the clinical responses of patients and their cog- co-deletion in glioma [21], MYC amplification and TP53 nate PDXs [32]. Although these data are more time- mutations [22], or activating alterations in KRAS and intensive and expensive to generate, it is still feasible to BRAF [5]. Indeed, although strong oncogenic KRAS and establish large in vivo screenings, covering a wide diver- BRAF alterations are mutually exclusive in treatment- sity of tumor types and drugs. PDXs are thus a clinically naïve tumors, a deeper allele-specific analysis identified a relevant platform for pre-clinical pharmacogenomic significant co-occurrence of activating RAS alterations studies and represent a more accurate approach to and BRAFD594 mutations, among other co-occurring identify predictive biomarkers compared with the use of hotspot mutations within MAPK signaling genes (i.e., cancer cell lines [33]. Mateo et al. Genome Medicine (2020) 12:78 Page 3 of 23 Here, we present a computational strategy to uncover which had been treated with a median of 18 (IQR 14– and exploit driver alteration co-occurrence patterns in 20) drugs. The final dataset consisted on 3127 experi- PDXs. By comparing the molecular profiles of responder ments performed on 187 PDXs and 53 treatment and non-responder PDXs to a given drug, we identify responses, across 5 tumor types: BRCA (breast cancer, driver co-occurrence networks and use them as a new n = 38), CM (cutaneous melanoma, n = 32), COREAD type of drug-response indicator, applicable much beyond (colorectal carcinoma, n = 51), NSCLC (non-small cell known biomarkers. We apply our strategy to the largest lung carcinoma, n = 27), PAAD (pancreatic adenocarcin- panel of PDXs and drugs available to date [28] and pro- oma, n = 38), and 1 PDX without tumor type annotation. spectively validate our findings in vivo. Finally, we adapt our strategy to derive response predictive models dir- Molecular representativity of PDXs ectly from continuous clinical outcome measures, such We used the OncoGenomic Landscapes tool [36] to ob- as progression-free survival, and evaluate them on a co- tain a 2D representation of the molecular heterogeneity hort of breast cancer patients. of the 187 PDXs being analyzed, and compared it to that of large reference cohorts of cancer patients. We down- Methods loaded the precomputed 2D projections of the following Genomic data processing reference cohorts from the OncoGenomic Landscapes A total of 1075 PDX models were established as part of webserver (oglandscapes.irbbarcelona.org): PanCancer a large pharmacogenomics screening that used the “one (n = 15,212), BRCA (breast cancer, n = 2021), CM (cuta- animal per model per treatment” (1 × 1 × 1) experimental neous melanoma, n = 492), COREAD (colorectal carcin- design to assess the population responses to 62 treat- oma, n = 1442), LUAD (lung adenocarcinoma, n = 1486), ments [28]. We collected somatic mutations and copy LUSC (lung squamous cell carcinoma, n = 352), and number alterations for 375 of them and used the Cancer PAAD (pancreatic adenocarcinoma, n = 442). We Genome Interpreter resource [14] to classify protein- merged LUSC and LUAD samples in order to get a ref- coding somatic mutations and copy number variants erence cohort for NSCLC (non-small cell lung carcin- into predicted passenger or known/predicted oncogenic. oma) PDXs. We selected the 2D coordinates of the In order to increase the clinical translatability, we sub- subset of TCGA and MSKCC patients of each reference sampled both datasets to consider those oncogenic alter- cohort and represented their distribution in the PanCan- ations covered by MSK-IMPACT [34] or by Foundation cer landscape as a level plot using the 2D kernel density Medicine [35] targeted gene panels to obtain DCO net- estimate function of the “seaborn” python library with works and TCT4U models that could be directly used 20 levels and a color map that represents probability with those kind of molecular profiles, which are becom- density as heat in the background. We selected the 2D ing widely used in the clinical setting. coordinates of the 187 PDXs and represented their indi- vidual location with points, colored by tumor type. Drug-response data In the original dataset, a total of 62 treatment groups Drug-response prediction based on Cancer bioMarkers were tested in 277 PDXs across six indications. Drug re- database sponse was determined by analyzing the change in We manually mapped the set of 53 drugs and drug com- tumor volume with respect to the baseline along time. binations tested in the cohort of PDXs to the corre- They combined two metrics (Best Response and Best sponding drug families in the Cancer bioMarkers Average Response) into a modified RECIST classification database [14] using drug target information available in (mRECIST) with four classes: PD (progressive disease), ChEMBL and DrugBank (Additional file 1: Table S1). SD (stable disease), PR (partial response), and CR We successfully assigned 50 out of the 53 treatments, (complete response). For our analyses, we considered spanning 29 drug family annotations. We considered PDXs whose tumors progressed upon treatment (PD) as those genomic alterations showing a “complete match” non-responders and PDXs whose tumors stopped grow- with any of the reported predictive biomarkers and col- ing (SD) or regressed (PR, CR) as responders. After lapsed them at gene level. Note that we adopted a tissue applying this binary classification, we had to exclude 9 agnostic approach in the development of TCT4U, and treatments for which there were less than 5 PDXs in one therefore, we did not require that the tissue of origin of of the two response groups, lacking thus enough inter- the PDX matched the tissue or lineage in which each individual heterogeneity to model drug response. A total biomarker was identified. Nevertheless, this information of 276 PDXs were treated in at least one of the 53 treat- is provided in Additional file 1: Table S2 to enable other ment groups considered, each treatment being tested in researchers to perform stratified analyses. We consid- 29 to 246 animals, with a median of 43 (IQR 38–93). ered as “approved” biomarkers those ones that are cur- We could obtain the molecular profile for 187 of them, rently approved by the FDA or by the main clinical Mateo et al. Genome Medicine (2020) 12:78 Page 4 of 23 guidelines in the field, such as the National Comprehen- given treatment arm, we compared the observed co- sive Cancer Network (NCCN), the College of American occurrence rate to the random expectation under a null Pathologists (CAP), the Clinical Pharmacogenetics Im- model with preserved sample- and gene-wise alteration plementation Consortium (CPIC), or the European Leu- rates. We obtained this null model by generating 1000 kemiaNet guidelines. We considered the rest of random permutations of the genomic alteration matrix biomarkers, with varying supporting evidence, as “ex- with the R package BiRewire [38]. We computed the perimental” biomarkers. The Cancer bioMarkers data- average probability that the co-alteration rate observed base usually reports more than one biomarker per drug in the actual dataset is larger than the co-alteration rate or drug family, and often a single patient (or PDX) har- observed in the permuted datasets with methyl_diff [37]. bors several biomarkers of response and/or non- When the average probability was larger than 95%, we response for the same drug or drug family. We grouped considered that the pair of drivers showed a tendency response and non-response biomarkers at gene level and towards co-occurrence. calculated the balanced accuracy (BAcc; average between Additionally, we computed the probability of the dif- sensitivity and specificity) of the prediction made by ferential co-occurrence rate between responder and each gene in each treatment arm. non-responder PDXs. We computed the same probabil- We weighted the binary predictions made by each ity under the null model and compared its distribution gene and combined them to obtain a final prediction per to the probability of differential co-occurrence observed treatment and PDX (wCombbmk). in the actual dataset. We selected the following sets of pairs per treatment arm: (i) Resp_Ps are those pairs X X showing a significant tendency towards co-occurrence in wCombbmk ¼ i∈Resp BAcci ∙si − j∈NonResp BAcc j ∙s j responder PDXs and showing a 95% probability of being ð1Þ co-altered more often in responder than in non- responder PDXs; (ii) NonResp_Ps, which show a signifi- Resp: set of genes with or without predictive cant tendency towards co-occurrence in non-responder biomarkers of response (Si, binary) PDXs and a 95% probability of being co-altered more NonResp: set of genes with or without predictive often in non-responder than in responder PDXs; and biomarkers of non-response (Sj, binary) (iii) General_Ps, which show a significant tendency to- BAcc: balanced accuracy of the predictive biomarker in wards co-occurrence in the whole treatment arm. Add- a given treatment arm itionally, we required that the selected pairs were altered more than once in the corresponding group, with a Driver co-occurrence networks minimum inferred alteration rate of 5%. In the case of Differentially altered drivers Resp_Ps and NonResp_Ps, we additionally required that For each treatment, we aimed at identifying differentially the probability of differential co-occurrence rate was lar- altered driver genes (DiffD) across responder and non- ger than the 95% percentile of the distribution of prob- responder PDXs. To this end, we used methyl_diff [37], abilities obtained when comparing permuted samples. an analytical solution to estimate the probability of the inequality between the mutation rate of each driver gene Driver co-occurrence networks across response groups, modeled using beta distribu- The differentially altered drivers (General_DiffD, Resp_ tions. We identified three sets of genes per treatment DiffD, NonResp_DiffD) and pairs of co-altered drivers arm: (i) Resp_DiffD are those genes with more than 95% (General_Ps, Resp_Ps, and NonResp_Ps) can be probability of showing higher alteration rate in re- expressed in terms of co-occurrence networks, in which sponder PDXs than in non-responders, (ii) NonResp_ nodes representing differentially altered driver genes DiffD are those genes with more than 95% probability of (DiffD) or driver genes involved in a pair of co-altered showing higher alteration rate in the non-responder than drivers (DiP) are connected according to significant co- in responder PDXs, and (iii) General_DiffD are those occurrences (Ps). For each treatment arm, we obtained genes with more than 95% probability of showing differ- three of such networks: (i) a general network (General_ ential alteration rate between the two response groups. DCO), (ii) a responder network (Resp_DCO), and (iii) a Additionally, we required that the selected genes were non-responder network (NonResp_DCO). altered more than once in the corresponding group, with a minimum alteration rate of 5%. TCT4U drug-response classifiers We described the DCO networks with a matrix of Bool- Driver pairs ean vectors (1: altered, 0: unaltered) encoding the alter- To identify pairs of driver gene alterations occurring ation status of differentially altered drivers and drivers more often than expected in each response group of a participating in co-occurring pairs in each PDX (DiffD_ Mateo et al. Genome Medicine (2020) 12:78 Page 5 of 23 DiP). We put together all those vectors in the form of a Next, we assessed the specific contribution of co- matrix and used it to train a decision tree-based gradient occuring pairs of drivers to the prediction of drug re- boosting classifier (CatBoost [39]). We did not specify sponse. To this end, we ranked driver pairs by the the edges as features because CatBoost already considers strength of the interaction, which is natively computed first order interactions between all pairs of features, by CatBoost and used to generate drug-response predic- meaning that it natively exploits driver co-occurrences tions. For each pair of drivers A and B, we classified the to predict treatment outcome. We used 100 trees with a samples on the basis of the status of genes A and B. We maximum depth ranging from 1 to 7, a learning rate then computed the average SHAP value of samples ranging from 0.2 to 1, and a coefficient at the L2 within each of the four resulting categories. We repre- regularization term of the Logloss function ranging from sented the results as SHAP interaction plots to uncover 1 to 10. We chose the best set of hyperparameters based what is the effect of having a driver alteration in gene A on the on the AUC obtained in the fivefold cross- on the SHAP value of gene B, and vice versa. We per- validation of 30 iterations. Please, note that we repeated formed this analysis separately for pairs of drivers lo- the same procedure for each treatment arm with each of cated in the same chromosome and pairs of drivers that the three DCO networks described before (General_ are far apart in the genome. DiffD_DiP, Resp_DiffD_DiP, and NonResp_DiffD_DiP). We assessed the accuracy and robustness of each of the Experimental validation in PDXs three classifiers by performing an external leave-one-out We collected all the available molecular profiles of the cross-validation (LOOCV). We used the balanced accur- VHIO collection of breast cancer PDXs. Most PDXs acy of the LOOCV as weight to combine the General_ were profiled using a hybridization-based capture panel DiffD_DiP, Resp_DiffD_DiP, and NonResp_DiffD_DiP of 410 genes (MSK-IMPACT) [34]. As we did for the predictions generated for each drug-PDX instance into a training set, we used the Cancer Genome Interpreter re- final score, as described in Eq. 2. source [6] to filter out as many passenger alterations as possible. In the same way we did for the LOOCV, we de- scribed the molecular profile of each PDX according to   the DiffD_DiP feature vectors associated to each DCO wComb ¼ BAccCB General ∙ I fCBGeneral ¼1g  Pð2Þ s − I fCBGeneral ¼ − 1g  P r network and used them to predict the response to the þBAccCB Resp ∙I fCB Resp¼1g  P CB Resp − BAccCB NonResp ∙I 53 treatments in the TCT4U collection. For each PDX, fCB NonResp¼ − 1g  P CB NonResp we ranked all treatments based on the predicted re- sponse and focused on the predictions generated by the CBGeneral, CBResp, CBNonResp: CatBoost binary predic- 21 treatments that attained a balanced accuracy of at tions based on General_DCOs, Resp_DCOs, and Non- least 0.6 in the leave-one-out cross-validation. In order Resp_DCOs. Values of 1 and − 1 indicate response and to increase the novelty of our findings, we prioritized non-response. those predictions that were not in agreement with pre- PCB_General, PCB_Resp, PCB_NonResp: probability estimates dictions made by known predictive biomarkers. Of the for the predictions. drugs and PDXs available in our laboratories, we se- I{CB=j}: indicator function that takes a value of 1 when lected 8 positive and 6 negative predictions spanning the CatBoost predicts class j or a value of 0 otherwise. following treatments: MEK inhibitor (n = 2), Pi3K inhibi- BAccCB: balanced accuracy attained by a given CB clas- tor (n = 5), taxane (n = 2), Pi3K inhibitor + CDK4/6 in- sifier in the LOOCV. hibitor (n = 3), and CDK4/6 inhibitor (n = 2). We used exactly the same pipeline to obtain DCO- For each drug-PDX pair, 2 to 10 tumors were subcuta- based predictions, biomarker-based predictions, and pre- neously implanted in immunocompromised mice and dictions based on an integrated model (DCO + grown until they reached a volume of 120–150 mm3. biomarkers). Tumors were treated with either vehicle or the corre- sponding drug or combination at a clinically relevant Model interpretation dose. Tumor growth was measured at least twice per We computed the average feature importance in the week for approximately 20 to 40 days, when typically LOOCV and excluded from the DCO networks those tumor volume in the control group had doubled twice driver genes that did not contribute to the predictions. or more. Caliper measurements were converted into This effectively removed also a lot of pairs of co- tumor volume estimates using the formula (l · w · w) · (π/ occurring drivers that did not contribute either. We then 6), where l and w are the major and minor tumor axes, quantified the feature’s local contribution to each indi- respectively. The response was determined following the vidual prediction by analyzing the Shapley additive ex- mRECIST guidelines that were used in the PDX screen- planation (SHAP) values calculated within CatBoost. ing that we used as training set [28]. Basically, we Mateo et al. Genome Medicine (2020) 12:78 Page 6 of 23 calculated the percentage change in tumor volume from than 95% probability of showing higher alteration rate in baseline (ΔVolt = (Vt − Vi)/Vi · 100) and determined the the one third of patients showing the most durable or BestResponse as the minimum value of ΔVolt after 10 or shortest clinical benefit, respectively, compared to the more days of treatment. In order to capture tumor third of patients at the other extreme of the distribution. growth dynamics, we also calculated the BestAverageRe- Additionally, we selected as General_DiffD all those P sponse as the minimum value of 1=n  ni¼1 ΔVoli after genes with more than 95% probability of showing differ- 10 or more days of treatment. PDXs were classified into ential alteration rate between the two extreme popula- response groups according to the mRECIST criteria ap- tions. The same strategy was applied in the identification plied in the following order: of pairs of driver gene alterations occurring more often CR: BestResp < − 95% and BestAvgResp < − 40% than expected considering all patients (General_Ps) or PR: BestResp < − 50% and BestAvgResp < − 20% separately for the one third of patients that relapsed the SD: BestResp < 35% and BestAvgResp < 30% latest (Resp_Ps) or the earliest (NonResp_Ps). The PD: BestResp ≥ 35% and BestAvgResp ≥ 30%. remaining steps were applied exactly as described for the We used the R package Xeva [40] to generate the binary TCT4U methodology. In this setting with only tumor response plots shown in Fig. 5. one treatment per patient, high-confidence predictions were selected by optimizing the threshold of the global score to get a maximum false discovery rate of 30% in Adaptation of TCT4U to use continuous clinical outcome the LOOCV, which happened when we kept predictions measurements of response with a score above 0.23 and predictions of We obtained both genomic and clinical data for a total non-response with a score below − 0.26. of 216 patients with HR+/HER2− metastatic breast can- cer that were treated with a CDK4/6 inhibitor in com- Results bination with an aromatase inhibitor in metastatic Driver co-occurrence networks of drug response setting [41]. All patients underwent prospective clinical Although thousands of genomic profiles of patient tu- genomic profiling consisting on the identification of sin- mors are available, accurate information about pharma- gle nucleotide variants, small indels, and copy number cological interventions and treatment outcome has not alterations detected from matched tumor-normal se- been systematically collected [23] or has not been dis- quence data using the MSK-IMPACT targeted gene closed yet [24]. Thus, to bypass these limitations, we panels. We used the Cancer Genome Interpreter [14] to compiled drug-response data obtained in PDXs, since filter out passenger mutations and CNVs and keep only they preserve the overall molecular profile of the original known or predicted driver mutations or copy number al- tumor, and maintain its cellular and histological struc- terations. Detailed treatment history data was collected ture [32]. In particular, we based our study on 375 PDXs for each patient and included all lines of systemic ther- for which somatic mutations and copy number alter- apy from the time of diagnosis of invasive carcinoma to ations have been acquired, together with their response the study data lock in September 2017. The exact regi- to 62 treatments across six indications, using the “one men, as well as the dates of start and stop of therapy, animal per model per treatment (1 × 1 × 1)” experimental was also recorded. For the current analysis, we consid- design [28]. As suggested by the authors, we adopted the ered the treatment duration time as a measure of clinical Modified Response Evaluation Criteria in Solid Tumors benefit derived by patients whose biopsies were collected (mRECIST) [28, 42] to assess the change in tumor vol- prior to or within the first 60 days of therapy initiation. ume in response to treatment. We considered “re- We used the TCT4U model of response to ribociclib sponders” those PDXs that showed a complete response to predicted response to CDK4/6 inhibition, as described (CR), partial response (PR), or stable disease (SD), and before. Due to the differences in clinical outcome mea- “non-responders” those with a progressive disease (PD) surements between the training and the clinical cohort, status. we decided to adapt the TCT4U methodology to use Of the 62 drugs and drug combinations tested, we se- continuous clinical outcome measurements as training lected 53 treatment arms that showed significant inter- set, instead of binary classification of drug response individual heterogeneity (i.e., a sufficient number of based on tumor growth. Our strategy consisted on com- “responder” and “non-responder” tumors) to model drug paring extreme populations both to derive the DCO net- response. In total, these data comprised 3127 experi- works and to train the classifier. We partitioned the ments performed on 187 PDXs [28] for which we had, at population into three equally sized sets and applied the least, 5 responder and 5 non-responder PDXs. First, we methodology described above. In this exercise, we set assessed whether this set of PDXs is representative of the cutoffs at 4.2 months and 9.7 months. We selected as the genomic diversity observed in human tumors by Resp_DiffD or NonResp_DiffD those genes with more comparing their alterations to the oncogenomic profiles Mateo et al. Genome Medicine (2020) 12:78 Page 7 of 23 extracted from 13,719 cancer patients [36]. We found only worked with driver somatic mutations and copy that the 187 PDXs considered broadly covered the whole number alterations. For each treatment, we grouped oncogenomic landscape represented by the full cohort responder and non-responder PDXs, irrespective of (“PanCancer” cohort in Fig. 1). When analyzing tumor the origin of their tumors. We first identified driver types individually, we observed that, while the muta- alterations that showed a differential mutation rate tional diversity of some of them is perfectly reflected in between responder and non-responder PDXs. Next, the PDX samples (e.g., colorectal and cutaneous melan- we identified pairs of driver alterations occurring oma tumors), the distribution of mutated genes showed more often than expected in each subpopulation clear differences in others (e.g., NSCLC). As expected, under a null model that preserves both gene-wise and we observe that PDXs sharing tissue of origin are more sample-wise alteration rates. Moreover, we identified similar between them than to other PDXs and, more im- pairs of driver alterations that show a differential co- portantly, that the same level of similarity is maintained alteration rate between responder and non-responder between PDXs and patient samples (Fig. 1). Overall, PDXs. With the overrepresented drivers (nodes) and there are PDXs representing the most populated areas of pairs of co-occurring drivers (edges) identified, we the PanCancer cohort, suggesting that the full collection built a driver co-occurrence (DCO) network that of PDXs may be used in downstream analyses. characterizes responder PDXs, another DCO network We used the Cancer Genome Interpreter [14] to fil- that characterizes non-responder PDXs, and a third ter out passenger mutations from PDX profiles, and general one consisting of all drivers and co- Fig. 1 Molecular representativity of PDXs. OncoGenomic Landscape 2D representations of the molecular heterogeneity of the 187 PDXs annotated with both drug-response data and oncogenic alterations, compared with that of their corresponding reference cohorts of cancer patients from TCGA and MSKCC. The points represent the location of each individual PDX, colored by tumor type. The distribution of the 187 PDXs can be compared to the distribution of patient samples, represented as density color-scale map in the background: PanCancer, PAAD (pancreatic adenocarcinoma), COREAD (colorectal carcinoma), CM (cutaneous melanoma), NSCLC (non-small cell lung cancer), and BRCA (breast cancer). The boxplots show the proximity (median Jaccard similarity coefficient) of PDXs to the 5% nearest neighbors in each comparison. On the left, we show the clustering of PDXs based on tissue of origin by comparing the proximity of PDXs of a given tumor type among themselves and to PDXs of other tumor types. On the right, we show the clustering of PDXs with patient samples of the same tumor type compared to patient samples of other tumor types. Stars denote the p value of a Wilcoxon rank-sum test (*< 0.05, **< 0.01, ***< 0.001, and ****p value < 0.0001) Mateo et al. Genome Medicine (2020) 12:78 Page 8 of 23 occurrences associated with both treatment responses be used to predict treatment outcome. For each drug, (Fig. 2a). we used the DCO networks to statistically classify PDXs DCO networks for each of the 53 drugs are detailed in as non-responder or responder. The goal of this exercise Additional file 1: Table S2 and can be visualized using is to identify, among the available treatments, the best Cytoscape [43] (Additional file 2; https://doi.org/10.6084/ possible option for each individual based on its oncoge- m9.figshare.12789068.v1). The total number of drivers nomic profile. We thus named the set of developed and driver co-occurrences captured in the DCO networks drug-response classifiers Targeted Cancer Therapy for varied substantially among treatments, ranging from 8 to You (TCT4U). 89 driver genes (median of 29 nodes, IQR 15–49) and 6 to In brief, for each treatment arm, we combined the 177 pairs of drivers (median of 26 edges, IQR 16–93) probabilities assigned by three gradient boosting classi- overrepresented in PDXs treated with abraxane + gemcita- fiers (CatBoost), trained with response, non-response, bine and alpelisib (BYL719), respectively. However, when and general DCO networks, into a single prediction considering individual animals, the number of altered score per drug-PDX pair (Fig. 3a). To increase the clin- drivers and pairs of drivers was small and remained quite ical translatability of our approach, we repeated the cal- stable across treatments, with a median of only 4 genes culations considering only those alterations detectable (IQR 3–5) and 2 driver co-occurrences (IQR 1–3) per by the Memorial Sloan Kettering-Integrated Mutation PDX (Fig. 2b, Additional file 3: Fig. S1). Profiling of Actionable Cancer Targets (MSK-IMPACT) We next sought to assess the novelty of our DCO net- [34, 44] and the Foundation Medicine (FM) gene panels works by quantifying their overlap with the set of anno- [35], which contain probes to detect 410 and 287 mu- tated response/non-response biomarkers for each tated genes, respectively, and are widely used in clinical treatment [14] that fulfill the eligibility criteria applied to settings. Finally, we assessed the effectiveness of TCT4U build our DCO networks. Figure 2c shows that, although by comparing its predictive power to that of FDA- there is some overlap, our approach vastly expands the approved and experimental biomarkers (see the set of genes to be considered in downstream treatment “Methods” section for details). prioritization applications. In some treatments, we ob- We collected the change in tumor volume and the serve a significant overlap between DCO networks and mRECIST classification for a total of 3127 experiments known biomarkers (i.e., LEE011 + ENCORAFENIB, with reported treatment outcome, comprising 187 PDXs BYL719 + ENCORAFENIB, or the FGFR inhibitor tested for response to 53 treatments. Figure 3b shows LLM871, Additional file 1: Table S3), but in most treat- the predictive performance of the models in a leave-one- ment arms, the counts are too low to attain statistical out cross-validation setting, whereby the genomic profile significance. However, when aggregated across treat- of PDXs is used to predict response to each treatment. ments, the overall overlap is significant, with 48 out of We observe that TCT4U models are applicable to all 359 possible drug-driver associations captured among drug-PDX pairs (3127), while alterations in approved the 1856 total drug-driver pairs covered in TCT4U, from and experimental biomarkers can only be found in about the universe of 30,687 eligible pairs (Fisher’s test OR half of them (1758). However, wherever applicable, both 2.45, p value 1.91·10−7, see Additional file 1: Table S3). methods attain a similar overall accuracy. As many treat- Although we do not expect to recapitulate all previously ment arms do not have a balanced number of reported drug-driver associations, finding some of them responders and non-responders, we quantified the bal- is a good sign of functional relatedness between DCO anced accuracy separately for each treatment and calcu- features and known mechanisms of action. lated their average (Fig. 3c, Additional file 1: Table S4). We are aware that without having performed a stratified TCT4U and known biomarkers attained an average analysis, we cannot rule out the possibility that tumor balanced accuracy of 0.58 and 0.52, respectively. The lineage might be a source of indirect associations between balanced accuracy was greater than 0.6 for 21 treat- some genomic features and response to treatment, and we ment arms, yielding a total of 1317 high-confidence might miss interesting context-specific biomarkers. predictions. Although this subset of predictions only However, we believe that, whenever identified, those bio- covered 42% of all drug-PDX pairs, their average bal- markers that are less context sensitive and that are com- anced accuracy improved to 0.66 (Fig. 3c, Additional mon across different tumor types would be of special file 1: Table S4). interest due to their wider applicability domain. In the case of known and experimental biomarkers, they could only predict with a balanced accuracy higher TCT4U: a collection of 53 drug-response classifiers for than 0.6 the response to 4 of those 21 treatments, genome-driven treatment prioritization namely encorafenib, encorafenib + binimetinib, encorafe- We then explored whether the sets of differentially co- nib + buparlisib (BKM120), and paclitaxel. As for FDA- altered genes in responder and non-responder PDXs can approved biomarkers, they predicted drug response in Mateo et al. Genome Medicine (2020) 12:78 Page 9 of 23 Fig. 2 Computational strategy and description of driver co-occurrence (DCO) networks. a We inferred DCO networks from the analysis of 3127 in vivo experiments that screened the efficacy of 53 treatments against a panel of 187 molecularly characterized PDXs of several tumor types. We first compared the patterns of oncogenic mutations and CNVs in responder and non-responder PDXs, regardless of the tissue of origin of the tumors. Next, we identified sets of driver genes showing differential alteration rates between responders and non-responders (DiffD), which are represented as red or blue nodes in DCO networks, respectively. Additionally, we identified pairs of genes whose alteration co-occurred more often than expected under a null model with preserved gene-wise and sample-wise mutational frequencies (Ps), and that did so more often in one of the two response groups. We represented each pair of co-altered drivers as two nodes connected by an edge. We derived a responder, non-responder, and general DCO network for each treatment. b Gray bars show the number of drivers and pairs of co-occurring drivers included in each DCO network derived from whole exome sequencing data. Red boxplots represent the number of drivers or driver co-occurrences identified in each individual PDX. c Red and blue boxes represent the overlap between DCO drivers and genes with annotated biomarkers of response or non-response, respectively. We show in light colors the number of drivers in the DCO networks that were not previously associated to drug response. Likewise, gray bars indicate the number of previously known biomarker genes that were not included in our DCO networks. In this analysis, we only considered as eligible biomarkers those that were identified in two or more PDXs, which is a requirement that any driver needs to satisfy in order to be incorporated to a DCO network. Stars denote the p value of a Fisher’s exact test (*< 0.05, **< 0.01, ***< 0.001, and ****p value < 0.0001) Mateo et al. Genome Medicine (2020) 12:78 Page 10 of 23 Fig. 3 Targeted Cancer Therapy for You (TCT4U), a collection of drug-response classifiers based on DCO networks. a Given a new tumor sample, we compare it to the patterns of driver alterations and co-alterations associated to response or lack of response to any of the treatments in TCT4U, and rank the drugs accordingly, predicting whether a drug will or not be effective. By quantifying each feature’s local contribution to individual predictions (SHAP values), we can know which are the precise molecular determinants used by the classifier and use this information for functional interpretation of the predictions. b Waterfall plot representing 3127 in vivo experiments, sorted left to right from the worst to the best response of 187 PDXs to 53 treatments. The heatmaps below show the predictions of TCT4U in a leave-one-out cross-validation setting, and the predictions made on the basis of approved or experimental biomarkers. Each heatmap has three rows, which correspond to the predictions obtained when examining the whole exome (WES), IMPACT410 (IM), or Foundation Medicine (FM) targeted gene panels. The number of predictions and the average balanced accuracy (BAcc) across treatments are annotated along the y-axis. The subset of high-confidence TCT4U predictions corresponds to the ones yielded by models that attained a BAcc of 0.6 or higher. The subset of high-confidence biomarkers corresponds to the subset of clinically approved ones. c Predictive performance of TCT4U-, biomarker-, and DCO + biomarker-based models, all trained with CatBoost within the TCT4U pipeline. d The precision of each set of predictions is illustrated by the red and blue sections of the stacked bar plots, which represent the proportion of correct response and non-response predictions. Analogously, incorrect predictions are represented in faint colors. Missing predictions (NA) are represented in white to offer a comparative overview of the recall. e Stacked bar plots representing the precision and recall of all TCT4U predictions and all previously known biomarkers covered by WES profiles split by treatment arm Mateo et al. Genome Medicine (2020) 12:78 Page 11 of 23 210 of the drug-PDX pairs, spanning 59 PDXs and 13 Interpretation of DCO-based predictions treatments. Due to sample imbalance, the balanced In addition to the prediction performance of each model, accuracy could only be calculated for a subset of it is key to evaluate their potential to uncover patterns treatments, making an average of 0.56 (Fig. 3d, Add- that can generate new hypothesis and propose novel bio- itional file 1: Table S4). It is also remarkable that, markers. Current tree-based explanation methods allow even if they consider a much lower number of genes, us to understand how the model uses input features to both MSK-IMPACT- and FM-derived models make predictions. Beyond assigning the importance of achieved comparable prediction accuracies (Fig. 3b, each feature to the global prediction, we can compute Additional file 1: Table S4). the Shapley additive explanation (SHAP) values to quan- The comparison of TCT4U to the simple combin- tify the local contribution of each feature to the individ- ation of biomarkers showed an overall improvement ual predictions (i.e., how important is each driver gene with respect to the current standard. However, given for predicting drug response in a given PDX). the way we modeled the treatment decision setting, not As expected, in targeted therapies that are currently having a given biomarker does not contribute to the approved in biomarker-specific indications, the SHAP predictions (i.e., not having NRAS alteration does not values of their intended targets show that they usually predict for response to BRAF inhibition). In order to contribute to DCO-based predictions. For instance, in maximize the information extracted from previously re- most PDXs, BRAF alteration predicts for response to ported biomarkers and perform a more controlled com- binimetinib (MEK inhibitor) used alone or in combin- parison, we decided to use exactly the same pipeline ation with encorafenib (BRAF inhibitor). On the other implemented for TCT4U but training the classifiers hand, KRAS and NRAS alterations, which are upstream with the few features representing known biomarkers. of BRAF, contribute positively when predicting response Reassuringly, the resulting predictions were strongly to binimetinib but negatively when predicting response correlated with the predictions based on the simple to its combination with encorafenib (Fig. 4a). Those are combination of biomarkers (Spearman’s rho 0.32, p the kind of predictions that are successfully achieved by value 1.2·10−100, Additional file 3: Fig. S2A). Known simple biomarkers (Fig. 3c, Additional file 3: Fig. S2). biomarkers were only able to predict response to 6 However, it becomes clear that, even when looking at an treatments with a balanced accuracy greater than 0.60, FDA-approved biomarker such as BRAF, single driver even when using gradient boosting classifiers instead of gene alterations do not contribute equally in all PDXs their simple combination. It is noteworthy that the six because of the effect of their interaction with additional treatments consist on encorafenib used alone or in driver alterations. combination with alpelisib, buparlisib, binimetinib, and PIK3CA-mutant tumors are sensitive to isoform- cetuximab. Figure 3c shows that known biomarkers selective PI3K inhibitors such as alpelisib (BYL719) [45– outperform TCT4U in the only four treatments involv- 47]. Accordingly, we observed a higher response rate ing encorabenib (256 predictions in total), while (65%, 15 out of 23) among PDXs with oncogenic biomarker-based models did not perform much better PIK3CA alterations compared to PDXs with wild-type than random in the remaining 17 treatments (Add- PIK3CA (44%, 52 out of 117), and the average SHAP itional file 1: Table S4). value of PIK3CA in altered PDXs is positive (1.16, Add- Further, we built integrated models containing both itional file 1: Table S2). However, PIK3CA-independent DCO-based features and known biomarkers, which mechanisms of PI3K activation (e.g., activating alter- slightly outperformed both DCO-only and biomarkers- ations in PIK3CB or PTEN loss) often limit the response only models, suggesting that the two sets of features to this treatment [48, 49]. The alpelisib DCO network carry orthogonal information (Fig. 3c). However, the contains three proteins involved in PI3K signaling, overall performance taking into account all treatment namely PIK3CA, PIK3R1, and PIK3C2B. Interestingly, arms did not change much with respect to the DCO- we found that PIK3CA-altered PDXs having no co- based models (Additional file 1: Table S4, Additional file occurring oncogenic alterations in the PI3K pathway 3: Fig. S2B-C). showed an even higher response rate (79%, 11 out of 14) Finally, while the coverage of approved or experimen- than those with co-occurring alterations in PIK3R1 or tal biomarkers is mostly limited to BRAF/MEK inhibi- PIK3C2B (44%, 4 out of 9). This is reflected in the nega- tors, PI3K/mTOR inhibitors, or cell cycle-related tive average SHAP values observed for PIK3R1 and treatments, the predictions made by TCT4U also cov- PIK3C2B in altered PDXs (− 0.83 and − 0.12, Additional ered other drug families including chemotherapies, RTK file 1: Table S2). Our DCO-based model was able to cap- inhibitors, and more experimental treatments targeting ture alterations that are likely to activate PI3K signaling Wnt signaling (WNT974), or apoptosis-related pathways in a PIK3CA-independent manner, which could limit the (TAS266, LGW813) (Fig. 3e). response to alpelisib [49]. Mateo et al. Genome Medicine (2020) 12:78 Page 12 of 23 Fig. 4 (See legend on next page.) Mateo et al. Genome Medicine (2020) 12:78 Page 13 of 23 (See figure on previous page.) Fig. 4 Specific contribution of driver co-occurrences to the prediction of drug response. a Summary plot of the local feature’s contribution (SHAP values) attributed to individual genes when predicting response to treatment with TCT4U. Each point represents the contribution of a driver gene to the prediction of response in a given PDX. The color of the points indicates whether the given driver gene was altered or not in each PDX. The three examples represent the five most explanatory genes when predicting response to three approved targeted therapies with biomarker- specific indications. b SHAP interaction plots and driver co-occurrence (DCO) networks representing the oncogenic alterations and pairs of alterations that are overrepresented in responder and non-responder PDXs. SHAP interaction plots show the effect of having driver alteration in gene A on the distribution of SHAP values of gene B. Each point represents the average SHAP value of PDXs classified on the basis of the status of the two drivers that tend to be co-altered. The size of the points is proportional to the number of PDXs that belong to each of the four resulting categories. The figure represents feature interactions involving driver genes that are located far apart in the genome. We also show some exemplary DCO networks with the driver co-occurrences represented in the accompanying SHAP interaction plots highlighted in yellow. The size of the nodes represents the average feature importance in the LOOCV, and their color represents the probability of being overrepresented in responder (red) or non-responder (blue) PDXs. Previously known biomarkers detected in the cohort are annotated with diamond shapes We extended this analysis to systematically assess the queried through cBioPortal [50, 51] (log2 OR 2.21 and > interaction between all co-occurring driver alterations in 3, respectively; both p values < 0.001). It has been sug- the 21 treatment arms that yielded high-confidence pre- gested that the amplitude of the regions affected by copy dictions. We identified a total of 253 unique feature in- number changes strongly determines patient prognosis teractions between co-occurring drivers (Additional file [52, 53]. Broader amplifications are likely to modify the 1: Table S2). Topological analysis of DCO networks re- dosage of multiple genes which, based on our observa- vealed that several large, strongly connected modules tions, could have a different impact on drug response were composed of genes that tend to be co-amplified or than the focal alteration of a single driver gene. co-deleted as part of the same genomic segment. Al- However, when co-altered drivers are genomically though in most of the cases those clusters of co-altered linked, it is very difficult to disentangle which is the spe- drivers are probably hitchhiking with a major driver al- cific contribution of each alteration because they would teration that is positively selected, in some other cases, still co-occur even if only one of them was actually con- the co-amplification and simultaneous overexpression of tributing to differential drug response. For this reason, adjacent oncogenes can provide a basis for cellular we decided to distinguish between interactions involving cross-talk. genes located in the same chromosomal arm (Additional For example, the co-amplification of the two FGFR li- file 3: Fig. S3) from those involving genes that are far gands FGF3 and FGF4 (chr11:q13) contributes positively apart in the genome (Fig. 4b). to the prediction of response to alpelisib (average SHAP As mentioned above, BRAF-driven tumors show a 1.59), while the alteration of FGF4 alone contributes higher response rate (83%, 15 out of 18) when treated negatively (average SHAP − 1.34). Similarly, the co- with encorafenib + buparlisib (BKM120) than BRAF alteration of MET (chr7:q31) and its downstream signal wild-type tumors (36%, 5 out of 14), with an average transducer BRAF (in chr7:q34, 24 Mb away) contributes SHAP value of 1.43 in altered vs. − 2.06 in non-altered more negatively to predict response to the inhibition of PDXs. However, the SHAP value is even more positive Wnt signaling by WNT974 than any of the two alter- for the subset of tumors with co-alteration of PTEN and ations alone (Additional file 3: Fig. S3B, Additional file 1: BRAF (1.89). Indeed all tumors with this co-alteration Table S2). In the same direction, MET alteration has responded to the treatment. On the contrary, PTEN also a negative impact on the prediction of response to deficiency in BRAF wild-type tumors is a negative pre- encorafenib, encorafenib + buparlisib (BKM120), and dictor of response (average SHAP value of − 2.60). Al- encorafenib + binimetinib in BRAF-altered PDXs (Add- though to a lesser extent, we observed the same trend in itional file 3: Fig. S3A, Additional file 1: Table S2). On the encorafenib + binimetinib model (Fig. 4b). Import- the other hand, the same co-alteration contributes more antly, PTEN and BRAF are also co-altered more often positively to predict response to encorafenib and encora- than expected in cancer patients queried through cBio- fenib + binimetinib than BRAF or MET alterations alone Portal (OR 1.86, p value < 0.001). (Additional file 3: Fig. S3A, Additional file 1: Table S2). Another remarkable example is the interaction be- Interestingly, a feedback loop between those two genes tween HERC2 and KRAS in the prediction of response has been shown to influence response to BRAF and/or to BINIMETINIB-3.5MPK (MEKi). This regimen was MET inhibitors [45]. tested on 25 pancreatic tumors, most of which were Not surprisingly, the two aforementioned pairs of KRAS driven (24 out of 25). Despite this fact, the down- drivers are also co-altered more often than expected in a stream inhibition of MEK1/2 was only effective in 62.5% non-redundant set of 46,697 cancer patient samples of them (15 out of 24). Interestingly, HERC2 was co- Mateo et al. Genome Medicine (2020) 12:78 Page 14 of 23 altered with KRAS in a substantial fraction of non- We treated five PDXs with alpelisib, three of which responder KRAS-driven tumors (25%, 5 out of 20). (PDX131, PDX293, and PDX156) were predicted as re- Accordingly, HERC2-KRAS co-alteration contributes sponders to the drug by TCT4U models, and two as negatively to the prediction of response (average SHAP non-responders (PDX191 and PDX153). The three PDXs value − 3.22), whereas having wild-type HERC2 and al- predicted as responders showed a co-amplification of tered KRAS contributes positively (0.54). In patients, the FGF3-FGF4-CCND1 triplet, located in the 11q13 HERC2 and KRAS are also significantly co-altered (cBio- genomic segment. In our DCO models, this triplet hap- Portal OR 1.77, p value < 0.001). pened more often in responder than in non-responder Overall, those examples illustrate how TCT4U is able PDXs, with an alteration rate of 7.46% and 1.37%, re- to detect and exploit interactions between driver alter- spectively. The three genes contributed positively to the ations to predict drug response in PDXs. Sometimes, prediction of response in these PDXs, with average these interactions involve previously known biomarkers, SHAP values of 0.72, 1.86, and 0.26 (Additional file 1: as illustrated by the co-alteration of BRAF and PTEN, Table S5). It is worth noting that our model, which was which would have antagonistic effects in predicting re- derived from 140 PDXs of different tumor types (i.e., 38 sponse to the combined inhibition of BRAF and PI3K BRCA, 42 COADREAD, 25 NSCLC, and 35 PDAC), did signaling. Moreover, many of these co-occurrence not show a significant tendency towards co-occurrence patterns are also found in patient cohorts, indicating a of PIK3CA and the 11q13 amplicon (OR 2.69, p value potential clinical translation of these findings. We pro- 0.26). Dysregulation of FGFR signaling can lead to vide the full set of interactions and their impact on downstream activation of PI3K/AKT pathway, and in- SHAP values in Additional file 1: Table S2. deed, a recent study reported that 73% of patients (8 of 11) with both an alteration in the PI3K/AKT/mTOR pathway and FGF/FGFR amplification experienced clin- Experimental validation of TCT4U drug-response ical benefit when treated with therapy targeting the predictions on a prospective PDX dataset PI3K/AKT/mTOR pathway, whereas only 34% of pa- In addition to the in silico benchmarks, we sought to tients (12 of 35) with PI3K/AKT/mTOR alterations prospectively evaluate the performance of the TCT4U alone did so [54]. However, the implication of FGF sig- models in new tumors. To this aim, we selected, from naling with respect to the clinical benefit of PI3K/AKT/ our VHIO collection of molecularly characterized breast mTOR blockage remains controversial. The retrospect- cancer PDXs, a subset of 14 drug-PDX pairs, with 8 tu- ive analysis of a large subset of patients enrolled in the mors predicted to respond and 6 predicted not to re- BOLERO-2 trial [55] showed that alterations in FGF sig- spond, namely alpelisib, an isoform-selective PI3Kα naling had a negligible impact (FGFR1) or slightly de- inhibitor (BYL719, n = 5); ribociclib, a CDK4/6 inhibitor creased (FGFR2) the clinical benefit of everolimus (LEE011, n = 2); the combination of both (alpelisib + treatment. In line with these findings, ER+/ERBB2− ribociclib, n = 3); the MEK inhibitor binimetinib (n = 2); metastatic breast cancer patients with FGFR1 and and paclitaxel, a taxane (n = 2). This is a particularly FGFR2 amplification did not derive a clinical benefit challenging set of drug-PDX pairs since, except in one from alpelisib + letrozole [56]. Accumulating evidence case, the anticipated TCT4U outcome did not agree with suggests that FGF signaling induced by FGFR1/2 ampli- approved or experimental biomarkers, either because the fication attenuates the response to PI3K blockage in individual genomic profiles did not have any biomarker PIK3CA-mutant breast cancer. However, the impact of altered (n = 8), or because the TCT4U predictions were FGF signaling in response to alpelisib in PIK3CA wild- opposed to those suggested by known biomarkers (n = type tumors originated from breast, as well as from 5). other tissues, has yet to be determined. We subcutaneously implanted the tumors in immuno- In our dataset, the three PDXs responded to the treat- compromised mice and let the tumors grow until they ment. In particular, in PDX293, we observed a partial re- reached a volume of 120–150 mm3. We then treated the sponse (PR) after 18 days of treatment, with a reduction PDXs for 15–57 days and measured their response to of 65% in the initial tumor volume. PDX131 and the administered drugs following the mRECIST guide- PDX156 showed a stable disease (SD) after 20 and 11 lines (see the “Methods” section for details). The days of treatment, respectively. On the other hand, complete results of our study, including TCT4U predic- PDX191 was predicted to be non-responder because, in tions, known biomarkers, treatment setting (drug dose, addition to the aforementioned 11q13 amplicon, it had duration, etc.), and tumor response (tumor growth, driver alterations in FGFR1, MYC, and GNAS that were mRECIST classification, etc.) for every PDX, can be contributing negatively to the prediction of response found in Additional file 1: Table S5 and are summarized (SHAP values of − 1.81, − 1.04, and − 0.73). In agreement in Fig. 5. with our prediction, the tumor increased its volume by Mateo et al. Genome Medicine (2020) 12:78 Page 15 of 23 Fig. 5 Experimental validation of TCT4U predictions. a The waterfall plot summarizes the results of 14 in vivo experiments comprising 5 PDXs treated with the alpelisib (BYL719) isoform-selective PI3K inhibitor (PI3Ki), 2 PDXs treated with ribociclib (LEE011) CDK4/6 inhibitor (CDKi), 3 PDXs treated with the combination of both (PI3Ki + CDKi), 2 PDXs treated with binimetinib MEK inhibitor (MEKi), and 2 PDXs treated with paclitaxel. The heatmap below shows the predictions based on TCT4U compared to previously known biomarkers. b–e Exemplary tumor response curves and features in the DCO network explaining TCT4U predictions. The size of the nodes represents the average feature importance in the LOOCV. The color of the nodes represents the contribution of each feature to the prediction of response in the given tumor (SHAP values). Previously known biomarkers detected in a given tumor are annotated with diamond shapes. The complete results of our study are provided in Additional file 1: Table S5 80% after 13 days of treatment (PD) (Additional file 1: alterations that together showed a stronger statistical as- Table S5). PDX153 was the only PDX with an oncogenic sociation than PIK3CA status, although they proved to PIK3CA mutation (p.K111E) reported to confer sensitiv- be less informative. ity to the treatment [14], and indeed, we observed a We administered ribociclib, a CDK4/6 inhibitor, to significant reduction of 83% in the tumor volume after PDX4 and PDX244_LR1, with the TCT4U prediction 35 days of treatment (i.e., a PR outcome). Our model that the two tumors would not respond to the drug. classified this PDX as non-responder because it also had PDX4 did not present any known biomarker of drug re- other alterations overrepresented among non-responder sponse, and the most influential feature in the ribociclib PDXs, such as MAP2K4 (SHAP value − 2.85) or NCOR1 DCO model was GNAS amplification. While the co- (SHAP value − 1.68). The DCO networks also consid- alteration of GNAS and AURKA is usually predictive of ered PIK3CA status, which tends to be more frequently response to ribociclib (average SHAP value 0.461), the altered in responder PDXs (22.29%) than in non- alteration of GNAS in the absence of additional alter- responder PDXs (11.40%; SHAP value 1.86). However, ations in this specific PDX contributed negatively to the the final prediction was driven by additional oncogenic prediction, with a SHAP value of − 0.88 (Additional file Mateo et al. Genome Medicine (2020) 12:78 Page 16 of 23 1: Table S5). On the other hand, we also treated contributed to the prediction of non-response in PDX244_LR1, which is a model of acquired resistance to PDX270, with SHAP values of − 2.147 and − 0.330. MYC ribociclib derived from a responder parental tumor was also altered in PDX288, but in this case, it was co- (PDX244). Accordingly, PDX244_LR1 simultaneously altered with SOX17. This co-alteration is distinctive of showed known biomarkers of response (CDKN2A- responder PDXs, with an observed co-occurrence rate of CDKN2B co-deletion) and non-response (TP53 13.76% in responders and 7.55% in non-responders. p.C176R) to the treatment [14]. TP53 was not included FGFR1 alteration also contributed to the prediction of in the ribociclib DCO network because we did not find response. In this case, neither tumor presented known it to be differentially altered between responders and biomarkers of response to MEK inhibition. When non-responders (44.65% vs. 48.39%). In line with what treated with binimetinib, PDX270 was classified as non- has been reported, CDKN2A-CDKN2B co-deletion was responder (PD), as the tumor volume had increased by slightly more common in responder than in non- 144%, even more than in untreated animals (117%). On responder PDXs (32.78% vs. 25.34) and contributed the contrary, and validating the TCT4U models, positively to the prediction with SHAP values of 0.116 PDX288 responded well to treatment (SD), and tumors and 0.249 (Additional file 1: Table S5). However, did not show any significant growth. Interestingly, an in- PDX244_LR1 presents an oncogenic mutation in RB1 tegrative genomics screen performed in 229 primary in- (p.M695Nfs*26), which showed a strong association with vasive breast carcinomas identified the co-amplification lack of response to CDK4/6 inhibition in the DCO net- of MYC and the 8p11-12 cytogenetic bands, together works (4.29% vs. 12.65%) and is the most negative ex- with aberrant methylation and expression of several plainer in this model, with a SHAP value of − 1.226. RB1 genes spanning the 8q12.1-q24.22 genomic region [58]. is the primary target of CDK4/6, and its status is a key This observation coincides with our DCO network de- determinant of CDK4/6 inhibition efficacy [57]. Accord- rived from whole exome sequencing data, where we ingly, RB1 overexpression is reported to confer sensitiv- could detect the co-amplification of a large cluster of ity to CDK4/6 inhibition in prostate cancer, although its genes located in the 8p11-p12 (HOOK3, FGFR1) and loss or deletion is not currently reported as a non- 8q11.23-q24.22 genomic regions (TCEA1, SOX17, response biomarker [14]. Our experiments showed that, CHCHD7, NCOA2, COX6C, MYC, NDRG1) in and in agreement with TCT4U predictions, the tumors responder PDXs, but not in non-responder ones increased their volume between 45 and 215%, being thus (Additional file 1: Table S2, Additional file 3: Fig. S3). cataloged as PD. Finally, we explored the TCT4U prediction capacity in We also treated three PDXs (PDX173, PDX98, and cytotoxic chemotherapy, where specific oncogenic char- PDX39) with the same PI3Kα and CDK4/6 inhibitors in acteristics should be less related to treatment efficacy. combination (alpelisib + ribociclib). The three of them We selected PDX222 and PDX39 to be treated with pac- had oncogenic mutations in TP53 (p.R249S, p.R249S, litaxel. While PDX222 did not present any known bio- and p.V157I), which are associated with non-response to marker of response, PDX39 sowed an MCL1 CDK4/6 inhibition [14]. However, DCO networks found amplification, which has been reported to promote re- additional response-associated genomic features (i.e., sistance to anti-tubulin chemotherapeutics [14, 59]. Al- MYC in PDX98 with SHAP value of 0.558, or H3F3A in though PDX222 showed alterations that are slightly PDX39 with SHAP value of 0.321, Additional file 1: more common in non-responder than in responder Table S5), and thus, TCT4U models predicted them to PDXs (SOX17 and TP53, SHAP values − 0.177 and − respond to this drug combination. We found that, in- 0.208), it also presented an ERBB2 amplification that in deed, all three tumors responded to the combination our model contributed positively to the final prediction treatment: PDX173 became completely tumor free (CR), (SHAP value 2.885). Regarding PDX39, TP53 and PDX39 showed a reduction of 47% (SD), and PDX98 of H3F3A alterations were the main negative contributors 25% (SD). to the prediction of response. When treated with pacli- Two PDXs were treated with the MEK inhibitor bini- taxel, both PDXs showed a progressive disease (PD). metinib, with TCT4U models predicting PDX270 to be Overall, TCT4U models correctly predicted the out- non-responder and PDX288 to respond to the drug. come of 12 of the 14 treatments tested whereas known Both PDXs presented RB1 loss (a loss-of-function muta- biomarkers only predicted correctly 2 of the 14 treat- tion p.Y321* and a deletion, respectively), which is over- ment outcomes. In particular, one of the TCT4U represented in the non-responder DCO network (5.50% misclassified responses was correctly predicted by vs. 17.46%) and contributed negatively to the prediction known biomarkers, while the rest were either incorrect (SHAP value of − 0.549 in both cases). Additional alter- (4 of 6) or missing (8). In PDXs treated with alpelisib, ations necessarily contributed to the divergent prediction we correctly classified 3 out of 4 responders and 1 out of of these PDXs. Alterations in CCND3 and MYC 1 non-responders, which makes a balanced accuracy of Mateo et al. Genome Medicine (2020) 12:78 Page 17 of 23 0.875. We could not estimate the average balanced ac- the treatment as non-responders, and the 40 for which curacy across treatments due to the limited sample size the time to progression was longer than 9.7 months as and the imbalanced number of responder and non- responders. The resulting DCO networks for this treat- responder PDXs per treatment arm. Therefore, the pre- ment, which are relatively small compared to TCT4U dictions presented and the associated genomic features DCO networks, contain a total of 7 drivers and 6 co- remain hypothesis in need for further validation. occurrences (see Fig. 6a). We then used the DCO networks to derive the corre- Bringing TCT4U from the workbench to the clinics sponding TCT4U models, which should be able to pre- To explore the clinical potential of TCT4U method- dict whether a given patient will obtain a significant ology, we analyzed a cohort of 216 metastatic breast can- clinical benefit. In a leave-one-out cross-validation, cer patients being treated at the Memorial Sloan TCT4U models yielded confident scores for 95 out of Kettering Cancer Center [41], and for which we have re- the 216 patients in the cohort (see the “Methods” sec- corded information of their oncogenomic profile and tion). Of these, we predicted that 50 patients would re- clinical outcome (Additional file 1: Table S6). These spond and 45 would not respond to the treatment. A metastatic patients had received between 1 and 17 Kaplan-Meier analysis of the cross-validation showed rounds of treatments (median of 2) before being selected that patients predicted to relapse early, with a median for a trial to test a combination of CDK4/6 and aroma- time to progression of 5.4 months, derived little clinical tase inhibitors. Each tumor was genetically profiled, benefit compared to the 50 patients predicted to relapse using the MSK-IMPACT panel, and the clinical outcome later, whose median time to progression was significantly of the treatment was recorded as progression-free sur- longer (13.5 months, log-rank test p value 0.002, Fig. 6b). vival (PFS). In this study, one third of the patients did We obtained consistent results when fitting a Cox pro- not derive a clinical benefit and relapsed before 5 portional hazards regression model (correlation coeffi- months. At the other extreme of the distribution, one cient − 0.61, p value 0.002), indicating that TCT4U third of the patients could be treated for more than 10 scores are negatively associated with risk of relapse. The months and were considered to present a durable clin- performance of TCT4U models clearly surpasses that of ical benefit. We are aware that a threshold of 10 months known biomarkers for this drug combination. Although might not be relevant in a first line treatment setting, 54% (116 of 216) of patients had at least one annotated where this drug combination has shown to achieve a biomarker, which is a good coverage compared to other median PFS of 24 months [60]. However, the PFS de- treatments, we could not find a significant association creases in subsequent lines of therapy, and in a meta- between observed and predicted outcomes, at least in static setting where over half of patients have received terms of PFS (Fig. 6b). prior therapies, a PFS of more than 10 months might As for the interpretation of the DCO-based model, still be a good surrogate measure of the clinical benefit. MYC and MAP3K1 alteration were the strongest ex- We did not have PDXs treated with a combination of plainers of TCT4U predictions of response, with average CDK4/6 and aromatase inhibitors, and the closest SAHP values of 3.94 and 2.84 in altered patients. Note- TCT4U model for it was derived in response to CDK4/6 worthy, MYC alteration was in common between this inhibition (ribociclib), based on 71 responder and 100 model and the ribociclib model derived from PDXs. On non-responder PDXs. Using this model, only 22.7% the other hand, ESR1 and CCND1 were the strongest patients (49 out of 216) were predicted to respond to negative contributors to the predictions, with average treatment and the remaining 77.3% were predicted as SHAP values − 1.46 and − 0.98 (Fig. 6c). Oncogenic mu- non-responders. The majority of patients (78%) relapsed tations in ESR1 are common in metastatic and pre- within the first year of treatment, but unfortunately, we treated breast cancer, emerging as a mechanism of have no data in this clinical series as to whether these acquired resistance to endocrine therapies that can ul- tumors regressed, at least initially. It thus seems that the timately result in a lack of response to the combinational outcome measure used to train the TCT4U model therapy [61]. (mRECIST), based on relative tumor growth, is not ap- Interestingly, the concomitant alteration of FGF3 and propriate in most clinical settings. FGF4 with CCND1 (11q13 amplicon) has a positive im- Without a model for this specific drug combination, pact on the SHAP value of CCND1-altered patients, and with the aforementioned differences in outcome which goes from − 1.13 to − 0.94 (Fig. 6d). This positive measures, we decided to adapt our methodology to clas- interaction is even stronger when the amplification sify patients based on the duration of the treatment be- spans PAK1 gene (in 11q14), reaching an average SHAP fore cancer relapsed. For this, we divided the cohort in value of − 0.75. Indeed, the FGF3-FGF4-CCND1 triplet three groups and considered the 40 patients for which tends to be significantly co-altered with PAK1 only in the tumors relapsed before 4.2 months after the start of patients that relapsed late (p values of 0.037, 0.038, and Mateo et al. Genome Medicine (2020) 12:78 Page 18 of 23 Fig. 6 Application of TCT4U to predict treatment outcome in a clinical cohort of HR+/HER2− metastatic breast cancer patients. a Driver co- occurrence (DCO) networks representing the oncogenic alterations and pairs of alterations that are overrepresented in patients that relapsed early (non-responders) or in patients that derived a durable clinical benefit (responders) from CDK4/6 inhibition combined with an aromatase inhibitor. The size of the nodes represents the average feature importance in the LOOCV. The color of the nodes represents the probability that alterations in a given gene are overrepresented in responders (red) or non-responders (blue). Previously known biomarkers detected in the cohort are annotated with diamond shapes. b The Kaplan-Meier analysis of progression-free survival (PFS). TCT4U high-confidence predictions are better able to discriminate between patients that would experience early and late relapse than known biomarkers, with a median time to progression of 5.4 and 13.5 months, respectively. c Summary plot of the SHAP values attributed to individual genes. Each point represents the contribution of a given feature to the prediction of response in a given patient. The color of the points indicates whether a given driver gene was altered or not in each patient. d SHAP interaction plot showing the positive effect of the co-alteration of FGF3, FGF4, and PAK1 with CCND1 (all in chr11q13-14). The size of the dots is proportional to the number of patients in each category 0.027, Additional file 1: Table S2). It is noteworthy that precision oncology. Indeed, current efforts now focus on CCND1 amplification is one of the biomarkers that has the analyses of oncogenomic patterns to identify action- been previously associated to response to CDK4/6 block- able alterations, drugs to modulate them, and bio- ade. Based on our observations, the genomic context of markers to monitor response. Of particular interest are CCND1 alteration seems to be relevant in relation to computational platforms such as OncoKB [13] or the CDK4/6 inhibition and deserves further investigation. Cancer Genome Interpreter [14], which not only identify Our results suggest that the proposed methodology oncogenic alterations and potential targets, but also esti- could be used to derive DCO networks and train pre- mate their potential clinical applicability. Most current dictive models from the kind of data obtained from in- strategies focus on the identification of a single vulner- terim analyses in oncological clinical trials. Moreover, ability (i.e., driver gene) whose activity can be modulated whenever the time to detect a clinical benefit is reason- by a drug. However, given the complexity and hetero- able, such the 10 months in this study, TCT4U models geneity in tumors, and the high connectivity between could be derived with the first patients and used in cellular processes, every cancer might respond differ- population enrichment strategies to establish the bases ently to a certain treatment, depending on its global for new recruitments in adaptive trials. oncogenomic profile. Indeed, the analysis of the mutational landscape of Discussion cancer has also uncovered the existence of mutual exclu- Cancer sequencing projects have unveiled hundreds of sivity and co-occurrence patterns among driver gene al- gene alterations driving tumorigenesis, enabling terations [16, 62]. Many computational tools have been Mateo et al. Genome Medicine (2020) 12:78 Page 19 of 23 developed to identify those combinatorial patterns ex- patterns from the data and to exploit them for drug- perimentally (i.e., via CRISPR-Cas9 screens [63, 64]) or response prediction. computationally [65–72]. Patterns of mutual exclusivity With the current implementation of TCT4U, we can arise from functional redundancy, context-specific present a collection of drug-response predictive models dependencies (i.e., tumor type- or subtype-specific driv- for 53 treatments belonging to 20 drug classes, including ing alterations), or synthetic lethality interactions. While targeted and more conventional chemotherapies. In a functional redundancy has been used to reveal unknown cross-validation setting, our drug-response models functional interactions [72], the synthetic lethality con- attained a global accuracy similar to that of approved cept has been very successfully applied to the identifica- biomarkers, but could be applied to twice as many sam- tion of novel therapeutic targets [63, 64] or rational drug ples, including drug classes for which no biomarker is combinations [64], and to the prediction of drug re- currently available. Moreover, in an in vivo prospective sponse in cell lines [64] and patients [71]. validation, our models correctly predicted 12 out of 14 Although less studied, driver co-occurrences are often responses to 5 drugs tested on 13 tumors. However, we interpreted as a sign of synergy, and in some cases, they could not estimate the average balanced accuracy across have shown to be functionally relevant [16–22]. How- treatments due to the limited number of PDXs per treat- ever, they have not yet been exploited for drug-response ment arm and the actual imbalance between responder prediction. With the methodology presented in this and non-responders. manuscript, we compared the mutational profiles of tu- Obviously, our approach also suffers from some limita- mors that responded or did not respond to a certain tions. Due to the lack of systematic reporting of treat- drug to define driver co-occurrence (DCO) networks, ment history of the patients enrolled in genomic studies which capture both genomic structure and putative [23], it is difficult to match response to a drug with indi- oncogenic synergy. We then used the DCO networks to vidual molecular profiles from clinical data. This practic- train classifiers to identify the best possible treatment ally impairs the systematic assessment of the prediction for each tumor based on its oncogenomic profile. accuracy in patients for computational frameworks like The development of tools for personalized treat- TCT4U, PanDrugs [73], PANOPLY [75], iCAGES [76], ment prioritization based on genomic profiles is an or other in silico drug prescription tools such as the active field of research. Recently, Al-Shahrour and Cancer Genome Interpreter [14] or OncoKB [13]. Ex- colleagues presented PanDrugs [73], an in silico drug perimental validation of computational approaches is prescription tool that uses genomic information, path- time-intensive and very expensive. Therefore, beyond way context, and pharmacological evidence to the thorough experimental validation presented in this prioritize the drug therapies that are most suitable for manuscript, only PanDrugs and PANOPLY predictions individual tumor profiles. PanDrugs goes beyond the were experimentally validated, although on a single case single-gene biomarker by taking into account the collect- study performed on a PDX model that was treated with ive gene impact and pathway context of the oncogenic al- 5 drugs (PanDrugs) or 2 drugs (PANOPLY). terations identified in a given patient. However, it Given the limited clinical representativity of drug combines clinical evidence with in vitro drug screening screens performed on cell lines [15, 27, 74], we relied on data gathered from cancer cell line panels, which have patient-derived xenografts (PDXs) to implement our limited clinical translatability [15, 27, 30, 74]. strategy and to identify biomarkers of drug response. Al- PANOPLY [75] is another computational framework though PDXs have shown a good level of agreement that uses machine learning and knowledge-driven net- with the course of disease evolution and treatment re- work analysis approaches to predict patient-specific sponse observed in patient tumors [30–33, 77, 78], they drug response from multiomics profiles. This tool present some important drawbacks, such as the eventual shows a great potential, but the method strongly de- loss of intra-tumoral heterogeneity [79, 80] or certain pends on whole genome and transcriptome patient engraftment bias [30, 81]. Additionally, we have to con- data, which is not routinely acquired in clinical prac- sider that PDXs might not completely recapitulate the tice. Other methods like iCAGES [76] have been de- influence of the tissue of origin in tumors that have been veloped mainly to identify patient-specific driver implanted subcutaneously in immunodeficient mice and genes from somatic mutation profiles, which are later whose stroma has possibly regressed and/or been re- used to prioritize drug treatments. However, iCAGES placed by mouse stroma, altering thus their subclonal only considers drugs that directly target the identified evolution and response to treatments [77, 82]. However, driver alterations based on current FDA prescription our strategy can be readily adapted to derive drug- guidelines. All those methods rely on prior know- response models from continuous clinical outcome mea- ledge, which is incomplete and biased, and have not sures, such as progression-free survival, which better been conceived to identify novel co-occurrence represent the data acquired during routine clinical Mateo et al. Genome Medicine (2020) 12:78 Page 20 of 23 practice and in clinical trials. Indeed, we derived re- a complementary source of information for clinical sponse models on a clinical cohort of breast cancer decision-making. metastatic patients being treated with a combination of CDK4/6 and aromatase inhibitors, showing a good cor- Conclusions relation with progression-free survival. We believe that the computational framework presented, Most importantly, TCT4U drug-response DCO which goes beyond the single-gene approach by exploit- networks are interpretable and provide clear hints to ing co-occurrence patterns, could represent a significant identify the potential mechanisms of response or non- advance towards the development of effective methods response present in each tumor. However, one key for personalized cancer treatment prioritization, with po- challenge in interpreting driver alteration co-occurrence tential applications in population enrichment strategies patterns is that they can also emerge without necessarily in the context of adaptive clinical trials. Overall, our being synergistic if a pair of genes is affected by a com- strategy represents an opportunity to accelerate the mon mutagenic process. This commonly happens when identification and validation of complex biomarkers with several oncogenes are co-amplified as part of the same the potential to increase the impact of genomic profiling genomic region, and for this reason, we have annotated in precision oncology. the candidate co-occurring drivers with the genomic lo- cation of its constituent genes and distance between them. Moreover, co-occurrence patterns can also emerge Supplementary information Supplementary information accompanies this paper at https://doi.org/10. as a result of the exposure to other mutagenic processes 1186/s13073-020-00774-x. that increase the mutational burden, the chromosomal instability, or that leave specific mutational signatures Additional file 1. This compressed folder contains all the [68, 69, 83]. Context- or tumor type-specific dependen- supplementary tables referenced throughout the manuscript, together with their description. cies can also be a source of indirect associations with Additional file 2. This compressed folder contains a cytoscape session drug response. Although those confounding factors can for each treatment with all the driver co-occurrence networks generated obscure the biological interpretation of the DCO net- from whole exome sequencing (WES), IMPACT410 (IM) and Foundation works, they certainly provide valuable information for Medicine (FM) gene panels. This information complements the informa- tion provided as Additional file 1: Table S2. This file is available at https:// drug-response prediction. Therefore, DCO networks are doi.org/10.6084/m9.figshare.12789068.v1 a valuable asset for hypothesis generation that need to Additional file 3. This pdf file contains the Supplementary Figures that be complemented with orthogonal sources of evidence, are referenced throughout the manuscript, together with their figure and functional validation will always be needed to dem- legends. Higher resolution figures and extended data are available upon reasonable request. onstrate synergy. We also showed that our methodology is well suited to work with any custom gene panel, provided that the Acknowledgements selected genes contribute to the differences in response The authors would like to thank IRB Barcelona and VHIO colleagues for providing feedback on the work, and Dr. David Torrents (ICREA-BSC) for to the drug being analyzed. As the cost of clinical mo- critically reading the manuscript. V.S. thanks Faye Su (Novartis Oncology) for lecular profiling continues to drop, it is very likely that providing study reagents alpelisib (BYL719) and ribociclib (LEE011). more types of data can be integratively analyzed to im- prove drug-response prediction. However, in order to Authors’ contributions L.M., M.D.-F., and P.A. conceived and designed the study. L.M. and M.D.-F ensure the clinical translatability of our method in the developed the methodology and applied it to analyze and interpret the short term, we decided to focus on well-supported data. A.G-O, M.P., J.A., M.B., and V.S. provided mice model data and helped oncogenic alterations that are readily detectable by cost- during its analysis and interpretation. M.S., P.R., and S.C. provided patients’ data and helped during its analysis and interpretation. L.M., M.D.-F., and P.A. effective methods in the clinical setting. We acknow- wrote the manuscript. V.S., P.R., and S.C. reviewed the manuscript. P.A. ledge that this is a very conservative decision, and we supervised the whole study. All authors read and approved the final accept that we might be missing biologically relevant in- manuscript. formation (i.e., non-coding alterations, methylation events, or expression changes). Indeed, current clinical Funding L.M. is a recipient of an FPI fellowship. P.A. acknowledges the support of the biomarkers for patient stratification are mostly based on Spanish Ministerio de Economía y Competitividad (BIO2016-77038-R), the the detection of histopathological, cytogenetic, and im- European Research Council (SysPharmAD: 614944), and the Generalitat de munohistochemical changes that are not always detect- Catalunya (VEIS 001-P-001). V.S. is a recipient of a Miguel Servet grant from ISCIII (CPII19/00033) and receives funds from AGAUR (2017 SGR 540). The able at DNA sequence level. For example, breast cancer PDX program is supported by a GHD-Pink (FERO foundation) grant to V.S., patient stratification strategies based on ER/PR and A.G.-O. and M.P. received a FI-AGAUR and a Juan de la Cierva (MJCI-2015- ERBB2 status have proven to be very informative, both 25412) fellowship, respectively. M.S., P.R., and S.C. acknowledge the support of the NIH grants P30 CA008748, RO1CA190642, and the Breast Cancer Re- in terms of prognosis and response to treatment [84]. search Foundation. Additionally, P.R. receives funds from the Breast Cancer Accordingly, TCT4U predictions should be regarded as Alliance. Mateo et al. Genome Medicine (2020) 12:78 Page 21 of 23 Availability of data and materials on oncology drug development: a meta-analysis of clinical trials leading to Data from the Novartis Institutes for BioMedical Research PDX encyclopedia FDA approval. J Natl Cancer Inst. 2015;107(11):djv253. (NIBR PDXE) is publicly available in the original publication [28] (DOI: https:// 4. Prasad V. Perspective: the precision-oncology illusion. Nature. 2016; doi.org/10.1038/nm.3954). Patients’ data is publicly available through the 537(7619):S63. cBioPortal for Cancer Genomics (http://www.cbioportal.org/study?id=breast_ 5. Chang MT, Bhattarai TS, Schram AM, Bielski CM, Donoghue MTA, Jonsson P, msk_2018) and is also deposited in the European Variation Archive (EVA). Chakravarty D, Phillips S, Kandoth C, Penson A, et al. Accelerating discovery The accession number for the sequencing data is PRJEB29597, https://www. of functional mutant alleles in cancer. Cancer Discov. 2018;8(2):174–83. ebi.ac.uk/eva/?eva-study=PRJEB29597. The accession number for the 6. Bailey MH, Tokheim C, Porta-Pardo E, Sengupta S, Bertrand D, deposited RNA sequencing data reported in this paper [41] is GEO: Weerasinghe A, Colaprico A, Wendl MC, Kim J, Reardon B, et al. GSE122088, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE122088. Comprehensive characterization of cancer driver genes and mutations. The remaining datasets supporting the conclusions of this article are Cell. 2018;173(2):371–85 e318. included within the article (and its additional files). 7. Rubio-Perez C, Tamborero D, Schroeder MP, Antolin AA, Deu-Pons J, Perez- Llamas C, Mestres J, Gonzalez-Perez A, Lopez-Bigas N. In silico prescription Ethics approval and consent to participate of anticancer drugs to cohorts of 28 tumor types reveals targeting Mice experiments were conducted following the European Union’s animal opportunities. Cancer Cell. 2015;27(3):382–96. care directive (86/609/CEE) and were approved by the Ethical Committee of 8. Senft D, Leiserson MDM, Ruppin E, Ronai ZA. Precision oncology: the road Animal Experimentation of the Vall d’Hebron Research Institute. Human ahead. Trends Mol Med. 2017;23(10):874–98. tumor samples for PDX were obtained after approval from the Ethics 9. Das S, Lo AW. Re-inventing drug development: a case study of the I-SPY 2 Committee of the Vall d’Hebron University Hospital. Written informed breast cancer clinical trials program. Contemp Clin Trials. 2017;62:168–74. consent was signed by all patients. The analysis of patient data was 10. Simon R. Critical review of umbrella, basket, and platform designs for approved by the Memorial Sloan Kettering Cancer Center Institutional oncology clinical trials. Clin Pharmacol Ther. 2017;102(6):934–41. Review Board (IRB), and all patients provided written informed consent for 11. Pallmann P, Bedding AW, Choodari-Oskooei B, Dimairo M, Flight L, tumor sequencing and review of patient medical records for detailed Hampson LV, Holmes J, Mander AP, Odondi L, Sydes MR, et al. Adaptive demographic, pathologic, and treatment information. The research designs in clinical trials: why use them, and how to run and report them. conformed to the principles of the Helsinki Declaration. BMC Med. 2018;16(1):29. 12. Thorlund K, Haggstrom J, Park JJ, Mills EJ. Key design considerations for Consent for publication adaptive clinical trials: a primer for clinicians. BMJ. 2018;360:k698. Not applicable. 13. Chakravarty D, Gao J, Phillips SM, Kundra R, Zhang H, Wang J, Rudolph JE, Yaeger R, Soumerai T, Nissan MH, et al. OncoKB: a precision oncology Competing interests knowledge base. JCO Precis Oncol. 2017;2017:10.1200/PO.17.00011. M.S. has received research funds from Puma Biotechnology, AstraZeneca, 14. Tamborero D, Rubio-Perez C, Deu-Pons J, Schroeder MP, Vivancos A, Rovira Daiichi-Sankio, Immunomedics, Targimmune, and Menarini Ricerche; is a co- A, Tusquets I, Albanell J, Rodon J, Tabernero J, et al. Cancer Genome founder of Medendi.org; and is on the advisory board of Menarini Ricerche. Interpreter annotates the biological and clinical relevance of tumor P.R. reports consulting/advisory board for Novartis and institutional research alterations. Genome Med. 2018;10(1):25. support from Illumina and GRAIL, Inc. S.C. has received research funds in the 15. Jaeger S, Duran-Frigola M, Aloy P. Drug sensitivity in cancer cell lines is not past from Novartis and Eli Lilly and ad hoc consulting honoraria from Novar- tissue-specific. Mol Cancer. 2015;14:40. tis, Sermonix, Context Therapeutics, and Revolution Medicines. The remaining 16. Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, Dimitriadoy authors declare that they have no competing interests. S, Liu DL, Kantheti HS, Saghafinia S, et al. Oncogenic signaling pathways in the cancer genome atlas. Cell. 2018;173(2):321–37 e310. Author details 17. Huun J, Lonning PE, Knappskog S. Effects of concomitant inactivation of 1 p53 and pRb on response to doxorubicin treatment in breast cancer cell Joint IRB-BSC-CRG Program in Computational Biology, Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and lines. Cell Death Discov. 2017;3:17026. Technology, Barcelona, Catalonia, Spain. 2Experimental Therapeutics Group, 18. Tu Q, Hao J, Zhou X, Yan L, Dai H, Sun B, Yang D, An S, Lv L, Jiao B, et al. Vall d’Hebron Institute of Oncology, Barcelona, Catalonia, Spain. 3Human CDKN2B deletion is essential for pancreatic cancer development instead of Oncology and Pathogenesis Program, Memorial Sloan Kettering Cancer unmeaningful co-deletion due to juxtaposition to CDKN2A. Oncogene. Center (MSKCC), New York, NY 10065, USA. 4Department of Pathology, MSKC 2018;37(1):128–38. C, New York, NY 10065, USA. 5Breast Medicine Service, Department of 19. Dembla V, Somaiah N, Barata P, Hess K, Fu S, Janku F, Karp DD, Naing A, Medicine, MSKCC and Weill-Cornell Medical College, New York, NY 10065, Piha-Paul SA, Subbiah V, et al. Prevalence of MDM2 amplification and USA. 6Growth Factors Laboratory, Vall d’Hebron Institute of Oncology, coalterations in 523 advanced cancer patients in the MD Anderson phase 1 Barcelona, Catalonia, Spain. 7Department of Biochemistry and Molecular clinic. Oncotarget. 2018;9(69):33232–43. Biology, Universitat Autònoma de Barcelona, Bellaterra, Catalonia, Spain. 20. Laroche-Clary A, Chaire V, Algeo MP, Derieppe MA, Loarer FL, Italiano A. 8 Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Combined targeting of MDM2 and CDK4 is synergistic in dedifferentiated Catalonia, Spain. 9CIBERONC, Barcelona, Spain. 10Breast Cancer Group, Vall liposarcomas. J Hematol Oncol. 2017;10(1):123. d’Hebron Institute of Oncology, Barcelona, Catalonia, Spain. 11Department of 21. Lauber C, Klink B, Seifert M. Comparative analysis of histologically classified Medical Oncology, Hospital Vall d’Hebron, Universitat Autònoma de oligodendrogliomas reveals characteristic molecular differences between Barcelona, Barcelona, Catalonia, Spain. subgroups. BMC Cancer. 2018;18(1):399. 22. Ulz P, Heitzer E, Speicher MR. Co-occurrence of MYC amplification and TP53 Received: 13 May 2020 Accepted: 11 August 2020 mutations in human cancer. Nat Genet. 2016;48(2):104–6. 23. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, Kovatich AJ, Benz CC, Levine DA, Lee AV, et al. An integrated TCGA pan- References cancer clinical data resource to drive high-quality survival outcome 1. Stockley TL, Oza AM, Berman HK, Leighl NB, Knox JJ, Shepherd FA, Chen EX, analytics. Cell. 2018;173(2):400–16 e411. Krzyzanowska MK, Dhani N, Joshua AM, et al. Molecular profiling of 24. Consortium APG. AACR project GENIE: powering precision medicine advanced solid tumors and patient outcomes with genotype-matched through an international consortium. Cancer Discov. 2017;7(8):818–31. clinical trials: the Princess Margaret IMPACT/COMPACT trial. Genome Med. 25. Guinney J, Saez-Rodriguez J. Alternative models for sharing confidential 2016;8(1):109. biomedical data. Nat Biotechnol. 2018;36(5):391–2. 2. Schwaederle M, Zhao M, Lee JJ, Eggermont AM, Schilsky RL, Mendelsohn J, 26. Iorio F, Knijnenburg TA, Vis DJ, Bignell GR, Menden MP, Schubert M, Aben Lazar V, Kurzrock R. Impact of precision medicine in diverse cancers: a meta- N, Goncalves E, Barthorpe S, Lightfoot H, et al. A landscape of analysis of phase II clinical trials. J Clin Oncol. 2015;33(32):3817–25. pharmacogenomic interactions in cancer. Cell. 2016;166(3):740–54. 3. Jardim DL, Schwaederle M, Wei C, Lee JJ, Hong DS, Eggermont AM, Schilsky 27. Gillet JP, Varma S, Gottesman MM. The clinical relevance of cancer cell lines. RL, Mendelsohn J, Lazar V, Kurzrock R. Impact of a biomarker-based strategy J Natl Cancer Inst. 2013;105(7):452–8. Mateo et al. Genome Medicine (2020) 12:78 Page 22 of 23 28. Gao H, Korn JM, Ferretti S, Monahan JE, Wang Y, Singh M, Zhang C, Schnell alpha-selective inhibition with alpelisib (BYL719) in PIK3CA-altered solid C, Yang G, Zhang Y, et al. High-throughput screening using patient-derived tumors: results from the first-in-human study. J Clin Oncol. 2018;36(13): tumor xenografts to predict clinical trial drug response. Nat Med. 2015; 1291–9. 21(11):1318–25. 47. Andre F, Ciruelos E, Rubovszky G, Campone M, Loibl S, Rugo HS, Iwata 29. Einarsdottir BO, Bagge RO, Bhadury J, Jespersen H, Mattsson J, Nilsson LM, H, Conte P, Mayer IA, Kaufman B, et al. Alpelisib for PIK3CA-mutated, Truve K, Lopez MD, Naredi P, Nilsson O, et al. Melanoma patient-derived hormone receptor-positive advanced breast cancer. N Engl J Med. 2019; xenografts accurately model the disease and develop fast enough to guide 380(20):1929–40. treatment decisions. Oncotarget. 2014;5(20):9609–18. 48. Nakanishi Y, Walter K, Spoerke JM, O'Brien C, Huw LY, Hampton GM, 30. Bruna A, Rueda OM, Greenwood W, Batra AS, Callari M, Batra RN, Pogrebniak Lackner MR. Activating mutations in PIK3CB confer resistance to PI3K K, Sandoval J, Cassidy JW, Tufegdzic-Vidakovic A, et al. A biobank of breast inhibition and define a novel oncogenic role for p110beta. Cancer Res. cancer explants with preserved intra-tumor heterogeneity to screen 2016;76(5):1193–203. anticancer compounds. Cell. 2016;167(1):260–74 e222. 49. Juric D, Castel P, Griffith M, Griffith OL, Won HH, Ellis H, Ebbesen SH, 31. Krepler C, Sproesser K, Brafford P, Beqiri M, Garman B, Xiao M, Shannan B, Ainscough BJ, Ramu A, Iyer G, et al. Convergent loss of PTEN leads to Watters A, Perego M, Zhang G, et al. A comprehensive patient-derived clinical resistance to a PI (3) Kalpha inhibitor. Nature. 2015;518(7538):240–4. xenograft collection representing the heterogeneity of melanoma. Cell Rep. 50. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, 2017;21(7):1953–67. Byrne CJ, Heuer ML, Larsson E, et al. The cBio cancer genomics portal: an 32. Pompili L, Porru M, Caruso C, Biroccio A, Leonetti C. Patient-derived open platform for exploring multidimensional cancer genomics data. xenografts: a relevant preclinical model for drug development. J Exp Clin Cancer Discov. 2012;2(5):401–4. Cancer Res. 2016;35(1):189. 51. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, 33. Byrne AT, Alferez DG, Amant F, Annibali D, Arribas J, Biankin AV, Bruna A, Jacobsen A, Sinha R, Larsson E, et al. Integrative analysis of complex Budinska E, Caldas C, Chang DK, et al. Interrogating open issues in cancer cancer genomics and clinical profiles using the cBioPortal. Sci Signal. precision medicine with patient-derived xenografts. Nat Rev Cancer. 2017; 2013;6(269):pl1. 17(4):254–68. 52. Solimini NL, Xu Q, Mermel CH, Liang AC, Schlabach MR, Luo J, Burrows AE, 34. Cheng DT, Mitchell TN, Zehir A, Shah RH, Benayed R, Syed A, Anselmo AN, Bredemeyer AL, Li MZ, et al. Recurrent hemizygous deletions in Chandramohan R, Liu ZY, Won HH, Scott SN, et al. Memorial Sloan cancers may optimize proliferative potential. Science. 2012;337(6090):104–9. Kettering-Integrated Mutation Profiling of Actionable Cancer Targets (MSK- 53. Smith JC, Sheltzer JM. Systematic identification of mutations and copy IMPACT): a hybridization capture-based next-generation sequencing clinical number alterations associated with cancer patient prognosis. Elife. 2018;7: assay for solid tumor molecular oncology. J Mol Diagn. 2015;17(3):251–64. e39217. 35. Frampton GM, Fichtenholtz A, Otto GA, Wang K, Downing SR, He J, Schnall- 54. Wheler JJ, Atkins JT, Janku F, Moulder SL, Stephens PJ, Yelensky R, Valero V, Levin M, White J, Sanford EM, An P, et al. Development and validation of a Miller V, Kurzrock R, Meric-Bernstam F. Presence of both alterations in FGFR/ clinical cancer genomic profiling test based on massively parallel DNA FGF and PI3K/AKT/mTOR confer improved outcomes for patients with sequencing. Nat Biotechnol. 2013;31(11):1023–31. metastatic breast cancer treated with PI3K/AKT/mTOR inhibitors. 36. Mateo L, Guitart-Pla O, Duran-Frigola M, Aloy P. Exploring the Oncoscience. 2016;3(5–6):164–72. OncoGenomic Landscape of cancer. Genome Med. 2018;10(1):61. 55. Hortobagyi GN, Chen D, Piccart M, Rugo HS, Burris HA 3rd, Pritchard KI, 37. Raineri E, Dabad M, Heath S. A note on exact differences between beta Campone M, Noguchi S, Perez AT, Deleu I, et al. Correlative analysis of distributions in genomic (methylation) studies. PLoS One. 2014;9(5):e97349. genetic alterations and everolimus benefit in hormone receptor-positive, 38. Gobbi A, Iorio F, Dawson KJ, Wedge DC, Tamborero D, Alexandrov LB, human epidermal growth factor receptor 2-negative advanced breast Lopez-Bigas N, Garnett MJ, Jurman G, Saez-Rodriguez J. Fast randomization cancer: results from BOLERO-2. J Clin Oncol. 2016;34(5):419–26. of large genomic datasets while preserving alteration counts. Bioinformatics. 56. Mayer IA, Abramson VG, Formisano L, Balko JM, Estrada MV, Sanders ME, 2014;30(17):i617–23. Juric D, Solit D, Berger MF, Won HH, et al. A phase Ib study of alpelisib 39. Dorogush AV, Ershov V, Gulin A. CatBoost: gradient boosting with (BYL719), a PI3Kalpha-specific inhibitor, with letrozole in ER+/HER2- categorical features support. arXiv:1810.11363. https://arxiv.org/abs/181 metastatic breast cancer. Clin Cancer Res. 2017;23(1):26–34. 0.11363. 57. Shapiro GI. Genomic biomarkers predicting response to selective CDK4/6 40. Mer AS, Ba-Alawi W, Smirnov P, Wang YX, Brew B, Ortmann J, Tsao MS, inhibition: progress in an elusive search. Cancer Cell. 2017;32(6):721–3. Cescon DW, Goldenberg A, Haibe-Kains B. Integrative pharmacogenomics 58. Parris TZ, Kovacs A, Hajizadeh S, Nemes S, Semaan M, Levin M, Karlsson P, analysis of patient-derived xenografts. Cancer Res. 2019;79(17):4539–50. Helou K. Frequent MYC coamplification and DNA hypomethylation of 41. Li Z, Razavi P, Li Q, Toy W, Liu B, Ping C, Hsieh W, Sanchez-Vega F, Brown multiple genes on 8q in 8p11-p12-amplified breast carcinomas. DN, Da Cruz Paula AF, et al. Loss of the FAT1 tumor suppressor promotes Oncogenesis. 2014;3:e95. resistance to CDK4/6 inhibitors via the hippo pathway. Cancer Cell. 2018; 59. Wertz IE, Kusam S, Lam C, Okamoto T, Sandoval W, Anderson DJ, Helgason 34(6):893–905 e898. E, Ernst JA, Eby M, Liu J, et al. Sensitivity to antitubulin chemotherapeutics is 42. Therasse P, Arbuck SG, Eisenhauer EA, Wanders J, Kaplan RS, Rubinstein L, regulated by MCL1 and FBW7. Nature. 2011;471(7336):110–4. Verweij J, Van Glabbeke M, van Oosterom AT, Christian MC, et al. New 60. Tanguy ML, Cabel L, Berger F, Pierga JY, Savignoni A, Bidard FC. Cdk4/6 guidelines to evaluate the response to treatment in solid tumors. European inhibitors and overall survival: power of first-line trials in metastatic breast Organization for Research and Treatment of Cancer, National Cancer cancer. NPJ Breast Cancer. 2018;4:14. Institute of the United States, National Cancer Institute of Canada. J Natl 61. Preusser M, De Mattos-Arruda L, Thill M, Criscitiello C, Bartsch R, Ruhstaller T, Cancer Inst. 2000;92(3):205–16. de Azambuja E, Zielinski CC. CDK4/6 inhibitors in the treatment of patients 43. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, with breast cancer: summary of a multidisciplinary round-table discussion. Schwikowski B, Ideker T. Cytoscape: a software environment for integrated ESMO Open. 2018;3(5):e000368. models of biomolecular interaction networks. Genome Res. 2003;13(11): 62. Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, Xie M, Zhang Q, 2498–504. McMichael JF, Wyczalkowski MA, et al. Mutational landscape and 44. Zehir A, Benayed R, Shah RH, Syed A, Middha S, Kim HR, Srinivasan P, Gao J, significance across 12 major cancer types. Nature. 2013;502(7471):333–9. Chakravarty D, Devlin SM, et al. Mutational landscape of metastatic cancer 63. Behan FM, Iorio F, Picco G, Goncalves E, Beaver CM, Migliardi G, Santos R, revealed from prospective clinical sequencing of 10,000 patients. Nat Med. Rao Y, Sassi F, Pinnelli M, et al. Prioritization of cancer therapeutic targets 2017;23(6):703–13. using CRISPR-Cas9 screens. Nature. 2019;568(7753):511–6. 45. Juric D, Janku F, Rodon J, Burris HA, Mayer IA, Schuler M, Seggewiss- 64. Szlachta K, Kuscu C, Tufan T, Adair SJ, Shang S, Michaels AD, Mullen MG, Bernhardt R, Gil-Martin M, Middleton MR, Baselga J, et al. Alpelisib plus Fischer NL, Yang J, Liu L, et al. CRISPR knockout screening identifies fulvestrant in PIK3CA-altered and PIK3CA-wild-type estrogen receptor- combinatorial drug targets in pancreatic cancer and models cellular drug positive advanced breast cancer: a phase 1b clinical trial. JAMA Oncol. 2019; response. Nat Commun. 2018;9(1):4275. 5(2):e184475. 65. Wu H, Gao L, Li F, Song F, Yang X, Kasabov N. Identifying overlapping 46. Juric D, Rodon J, Tabernero J, Janku F, Burris HA, Schellens JHM, Middleton mutated driver pathways by constructing gene networks in cancer. BMC MR, Berlin J, Schuler M, Gil-Martin M, et al. Phosphatidylinositol 3-kinase Bioinformatics. 2015;16(Suppl 5):S3. Mateo et al. Genome Medicine (2020) 12:78 Page 23 of 23 66. Szczurek E, Beerenwinkel N. Modeling mutual exclusivity of cancer mutations. PLoS Comput Biol. 2014;10(3):e1003503. 67. Kim YA, Madan S, Przytycka TM. WeSME: uncovering mutual exclusivity of cancer drivers and beyond. Bioinformatics. 2017;33(6):814–21. 68. Dao P, Kim YA, Wojtowicz D, Madan S, Sharan R, Przytycka TM. BeWith: a Between-Within method to discover relationships between cancer modules via integrated analysis of mutual exclusivity, co-occurrence and functional interactions. PLoS Comput Biol. 2017;13(10):e1005695. 69. Canisius S, Martens JW, Wessels LF. A novel independence test for somatic alterations in cancer shows that biology drives mutual exclusivity but chance explains most co-occurrence. Genome Biol. 2016;17(1):261. 70. Mina M, Raynaud F, Tavernari D, Battistello E, Sungalee S, Saghafinia S, Laessle T, Sanchez-Vega F, Schultz N, Oricchio E, et al. Conditional selection of genomic alterations dictates cancer evolution and oncogenic dependencies. Cancer Cell. 2017;32(2):155–68 e156. 71. Lee JS, Das A, Jerby-Arnon L, Arafeh R, Auslander N, Davidson M, McGarry L, James D, Amzallag A, Park SG, et al. Harnessing synthetic lethality to predict the response to cancer treatment. Nat Commun. 2018;9(1):2546. 72. Vandin F, Upfal E, Raphael BJ. De novo discovery of mutated driver pathways in cancer. Genome Res. 2012;22(2):375–85. 73. Pineiro-Yanez E, Reboiro-Jato M, Gomez-Lopez G, Perales-Paton J, Troule K, Rodriguez JM, Tejero H, Shimamura T, Lopez-Casas PP, Carretero J, et al. PanDrugs: a novel method to prioritize anticancer drug treatments according to individual genomic data. Genome Med. 2018;10(1):41. 74. Domcke S, Sinha R, Levine DA, Sander C, Schultz N. Evaluating cell lines as tumour models by comparison of genomic profiles. Nat Commun. 2013;4:2126. 75. Kalari KR, Sinnwell JP, Thompson KJ, Tang X, Carlson EE, Yu J, Vedell PT, Ingle JN, Weinshilboum RM, Boughey JC, et al. PANOPLY: omics-guided drug prioritization method tailored to an individual patient. JCO Clin Cancer Inform. 2018;2:1–11. 76. Dong C, Guo Y, Yang H, He Z, Liu X, Wang K. iCAGES: integrated CAncer GEnome Score for comprehensively prioritizing driver genes in personal cancer genomes. Genome Med. 2016;8(1):135. 77. Hidalgo M, Amant F, Biankin AV, Budinska E, Byrne AT, Caldas C, Clarke RB, de Jong S, Jonkers J, Maelandsmo GM, et al. Patient-derived xenograft models: an emerging platform for translational cancer research. Cancer Discov. 2014;4(9):998–1013. 78. Izumchenko E, Paz K, Ciznadija D, Sloma I, Katz A, Vasquez-Dunddel D, Ben- Zvi I, Stebbing J, McGuire W, Harris W, et al. Patient-derived xenografts effectively capture responses to oncology therapy in a heterogeneous cohort of patients with solid tumors. Ann Oncol. 2017;28(10):2595–605. 79. Villacorta-Martin C, Craig AJ, Villanueva A. Divergent evolutionary trajectories in transplanted tumor models. Nat Genet. 2017;49(11):1565–6. 80. Eirew P, Steif A, Khattra J, Ha G, Yap D, Farahani H, Gelmon K, Chia S, Mar C, Wan A, et al. Dynamics of genomic clones in breast cancer patient xenografts at single-cell resolution. Nature. 2015;518(7539):422–6. 81. Willyard C. The mice with human tumours: growing pains for a popular cancer model. Nature. 2018;560(7717):156–7. 82. Wang M, Yao LC, Cheng M, Cai D, Martinek J, Pan CX, Shi W, Ma AH, De Vere White RW, Airhart S, et al. Humanized mice in studying efficacy and mechanisms of PD-1-targeted cancer immunotherapy. FASEB J. 2018;32(3): 1537–49. 83. Alexandrov LB, Nik-Zainal S, Wedge DC, Aparicio SA, Behjati S, Biankin AV, Bignell GR, Bolli N, Borg A, Borresen-Dale AL, et al. Signatures of mutational processes in human cancer. Nature. 2013;500(7463):415–21. 84. Onitilo AA, Engel JM, Greenlee RT, Mukesh BN. Breast cancer subtypes based on ER/PR and Her2 expression: comparison of clinicopathologic features and survival. Clin Med Res. 2009;7(1–2):4–13. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.