ARTICLE https://doi.org/10.1038/s41467-022-30582-3 OPEN Mutational signatures are markers of drug sensitivity of cancer cells Jurica Levatić1, Marina Salvadores 1, Francisco Fuster-Tormo1,3 & Fran Supek 1,2 ✉ Genomic analyses have revealed mutational footprints associated with DNA maintenance gone awry, or with mutagen exposures. Because cancer therapeutics often target DNA 1234567890():,; synthesis or repair, we asked if mutational signatures make useful markers of drug sensitivity. We detect mutational signatures in cancer cell line exomes (where matched healthy tissues are not available) by adjusting for the confounding germline mutation spectra across ancestries. We identify robust associations between various mutational signatures and drug activity across cancer cell lines; these are as numerous as associations with established genetic markers such as driver gene alterations. Signatures of prior exposures to DNA damaging agents – including chemotherapy – tend to associate with drug resistance, while signatures of deficiencies in DNA repair tend to predict sensitivity towards particular ther- apeutics. Replication analyses across independent drug and CRISPR genetic screening data sets reveal hundreds of robust associations, which are provided as a resource for drug repurposing guided by mutational signature markers. 1 Genome Data Science, Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology, C/ Baldiri Reixac 10, 08028 Barcelona, Spain. 2 Catalan Institution for Research and Advanced Studies (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain. 3Present address: MDS Group, Josep Carreras Leukaemia Research Institute, Ctra de Can Ruti, Camí de les Escoles s/n, 08916 Badalona, Spain. ✉email: fran.supek@irbbarcelona.org NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 C ancer precision medicine draws on the presence of constitute markers of drug sensitivity. Quantifying somatic somatically acquired changes in the tumor, which serve as mutational signatures in cell line genomes is however difficult, predictive markers of response to drugs and other thera- because a matched normal tissue from the same individual is pies. Commonly these markers are individual genetic changes, typically not available and thus cannot be used to remove the such as driver mutations affecting oncogenes or tumor suppressor abundant germline variation. After filtering the known germline genes, or copy-number alterations thereof. Many commonly variants listed in population genomic databases31,32, somatic employed cancer drugs act by interfering with DNA synthesis or mutations are still greatly outnumbered by the residual germline maintenance or by damaging DNA. Therefore, the altered capa- variants (Fig. 1a, b), which may confound downstream analyses city of cancer cells to repair and/or replicate DNA is the basis of such as the inference of mutational signatures. We introduce a many classical therapies, such as platinum-based agents, and also method to infer somatic mutational spectra from cancer genomes recently introduced or upcoming therapies, such as PARP inhi- without a matched control sample, while adjusting for the resi- bitors or ATR inhibitors (reviewed in refs. 1–3). It is paramount to dual germline variation. We apply this to infer trinucleotide identify predictive markers that are associated with failures of mutation signatures in cancer cell line exomes, and identify DNA maintenance in cancer cells. associations with sensitivity to drugs and to genetic perturbation However, while DNA repair is often deficient in tumors, many across cell line panels. Replication analyses across independent DNA repair genes such as MLH1, MGMT, BRCA1, or ATM do not data sets indicated that mutational signatures are broadly commonly bear somatic mutations. Instead, they are commonly applicable markers of drug sensitivity, matching or exceeding inactivated epigenetically4–6, or by alterations in trans-acting common genomic markers such as oncogenic driver mutations or factors7, and so their deficiencies are difficult to predict from the copy number alterations. gene sequence. Additionally, germline cancer-predisposing variants commonly affect DNA repair genes8–10, however, pathogenicity of such variants is often challenging to predict. Because of the above, Results other types of molecular markers may be more useful to infer An ancestry-matching approach removes subpopulation- about failed DNA repair. This is exemplified in ”BRCAness” – a specific trinucleotide spectra to accurately infer mutation gene expression signature that suggests a deficient homologous signatures. A substantial amount of the germline variation in a recombination (HR) pathway, even in the absence of deleterious cell line exome cannot be removed by filtering based on minor genetic variants in the BRCA1/2 genes. variant frequency in population databases (Fig. 1b). Therefore we In addition to gene expression, mutational signatures–readouts devised an approach to measure the somatic trinucleotide of genome instability–can characterize DNA repair deficiencies. mutation spectrum – the input for the inference of mutational One common type of signature describes relative frequencies of signatures33 – while rigorously adjusting for the contamination somatic single-nucleotide variants (SNV) across different trinu- by the residual germline mutation spectrum. cleotide contexts. Certain mutational signatures were found to be Because mutational processes differ across human associated with failures in DNA mismatch repair (MMR) and HR populations34, there is potential for this to confound analyses of pathways11 as well as DNA polymerase proofreading12,13 and the somatic mutation spectrum. Given the high number of base excision repair (BER)14–16 and nucleotide excision repair residual germline variants post-filtering (Fig. 1b), even slight (NER)17 failures. Inducing DNA repair deficiencies in cancer cell differences in the germline spectrum can cause large deviations in lines is able to reproduce some of these signatures18–21. Other the observed spectrum, which is a mix of somatic and germline types of mutation signatures based on small insertions and variation. deletions (indels)9 and on structural variants22 are also starting to To address this, we implemented an ancestry-matching be introduced. procedure, looking up the individuals with a similar ancestry to Because mutational signatures describe the state of the DNA each cell line’s ancestry. In particular, we clustered the cell line repair machinery of a cancer cell, they may be able to serve as a exomes together with germline exome samples from the TCGA drug sensitivity marker. This is exemplified by a mutational sig- data set, grouping by principal components derived from nature associated with pathogenic variants in BRCA1 and BRCA2 common germline variation (Fig. 1e; Methods section). The genes11,23, thus identifying HR deficient tumors. The signature is TCGA individuals clustered with a cell line provided a base- common in ovarian and breast cancers, but genomic analyses line germline mutational spectrum, which can be subtracted from have detected it across other cancer types24,25, suggesting the the observed mutation spectrum to estimate the somatic mutation potential for broad use of drugs that target HR-deficient cells, spectrum. such as PARP inhibitors. To this end, genomics-based predictors We benchmarked our ancestry-matching procedure for that draw on mutational signatures of HR deficiency have been the accuracy of reconstructing the correct somatic mutation developed26,27. We propose that this principle may extend to spectrum in a cancer cell line exome. To this end, we used SNV other types of mutational processes, potentially revealing tumor calls from TCGA cancer exomes where the matched normal was vulnerabilities. ignored, thus simulating the mutation calls that would be Human cancer cell line panels provide an experimental model obtained from cell line genomes (see Methods section). We then for the diversity in tumor biology that is amenable to scaling-up. compared the reconstructed somatic SNV mutation spectrum Drug screens and genetic screens on large cell line panels28,29 to the true somatic spectrum, obtained by contrasting tumor have identified correlations between the sensitivity to a drug (or exomes with the matched healthy tissue exomes from the same to a genetic perturbation), and the genetic, epigenetic, or tran- individuals. scriptomic markers in the cell lines. Encouragingly, genetic Ancestry-matching improves over the commonly used strategy, markers known to have clinical utility (e.g. BRAF mutations that is simply filtering out known germline variants according to for vemurafenib, EGFR mutations for gefitinib, BCR-ABL population genomic databases (Fig. 1c, 1d and Supplementary fusion for imatinib sensitivity) are also evident in cell line panel Fig. 1d); error in somatic trinucleotide frequency spectrum data analyses30, suggesting potential for discovery of further (Methods section) is 68.8 versus 124.1 across all tissues, while for useful genomic markers. comparison, the error expected by ‘self-similarity’ via a bootstrap- Here, we used large-scale cell line data to investigate the resampling of mutations from the same tumor samples would be hypothesis that mutational signatures in cancer genomes 67.5, close to that obtained via ancestry-matching. 2 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE Fig. 1 Evaluation of the ancestry-matching method to infer somatic mutation spectra on exomes without a matched normal control. a, b Germline variants greatly outnumber somatic mutations in exomes of various tumor types (n = 52 BRCA, 33 KIRC, 53 GBM, 19 BLCA, 15 LUSC, and 67 LUAD cancer exomes) (a), also after attempting to filter out germline variants according to the minor allele frequency (MAF) of variants listed in the gnomAD database (n = 450 cancer exomes) (b). The center line of box plots denotes medians of data points and the box hinges correspond to the 1st and 3rd quartiles, while whiskers extend to 1.5× IQR from the hinges. Data points beyond the end of the whiskers are shown individually. c Error between the real somatic 96 tri- nucleotide profiles and the profiles obtained with the ancestry-matching procedure, after various numbers of clusters (based on principal components of common germline variants; see Methods) are considered (n = 450 cancer exomes). d Comparison of the ancestry-matching method (with the number of clusters set to 13), the baseline procedure (variant filtering by population MAF<0.001%), regressing out mutational signatures reported as related to germline variants (signatures 1 and 5, and SNP signature31,110), and the error expected by chance (estimated by bootstrapping mutations). P-values by two- sided Wilcoxon rank-sum test (n = 450 cancer exomes). e A schematic representation of the ‘ancestry matching’ procedure. For compactness, the X-axes on the mutation spectra illustrations list only a subset of mutation types. PCA, principal components analysis. Error bars in panels b–d are the standard error of the mean. Source data are provided as a Source Data file. We considered various numbers of population clusters for 13, versus 71.5 for 3 clusters; this improvement is modest and according to the error in reconstruction of the correct somatic so the 3-cluster solution may also provide a satisfactory baseline trinucleotide spectrum. Selecting three clusters, expectedly, for downstream mutational spectra analyses). recovers the major ethnicity groups (European, Asian and Encouragingly, comparing the 13 ancestry clusters sorted by African, Supplementary Fig. 1a) and further increasing the self-reported ethnicity (Supplementary Fig. 1b), intra-ethnicity number of clusters to 13 minimizes the error in reconstructing trinucleotide mutational profiles are more similar than the inter- true somatic trinucleotide mutation spectra (Fig. 1c; error is 68.8 ethnicity profiles (Supplementary Fig. 1c), and a PC analysis of NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 the trinucleotide spectra of the rare variants separates the major (Supplementary Fig. 6 shows examples). Finally, some SBS-CL may ethnicity groups (Supplementary Fig. 1d). This is consistent with reflect contamination with residual germline variation, as well as reports of differential mutagenic processes in the human germline with sequencing artifacts, similarly as was recently reported for across ancestral groups – for example, the HCC>HTC (H = not many SBS signatures recovered from tumor genomes9. G) variants were reported to be increased in Europeans, NCG>NTG mutations in Native Americans and NAC>NCN and TAT>TTT in some East Asians34–36. These reports, together Mutational signatures predict cell line drug response more with our benchmarking using simulations, support the use of accurately than oncogenic mutations or copy number altera- ancestry-specific baselines in inferring somatic mutational spectra tions. Genetic and epigenetic alterations in cancer cell lines are of unmatched cancer genomes, such as cell line genomes. often investigated as markers of sensitivity to chemical We applied the ancestry-matching methodology (Fig. 1e) to compounds29,30. We hypothesized that mutational signatures in a exome sequencing data of 1071 cancer cell lines37, yielding their cell line genome can serve as similarly informative markers of somatic trinucleotide spectra. On this data, we performed de novo drug sensitivity or resistance. We compared their predictive discovery using an NMF approach, broadly as described by ability to that of the markers commonly used to predict drug Alexandrov et al.33 (with certain modifications, see Methods response in cell lines: oncogenic mutations (in 470 cancer driver section), where we extracted those NMF solutions that resembled genes30), recurrent focal copy number alterations (CNAs at 425 previously reported tumor mutational signatures9 of single base genes30), and DNA methylation data at informative CpG islands substitutions (SBS). We tested a number of variations on the data (HypMet at 378 genes30). Additionally, we examined gene filtering and the mutation extraction methodology (Supplemen- expression patterns (mRNA levels of 1564 genes that are either tary Data S1) to improve agreement with the known set of SBS represented in the L1000 assay39 or are known drug target signatures9 and their known distribution across tissues, as well to genes40), because gene expression can be highly predictive of drug improve power of the set of mutational signatures to predict drug response30,41, possibly because it reflects differences between responses in the cell lines (Methods section; Supplementary various cancer types and subtypes. Data S1). To further demonstrate the utility of the ancestry We predicted the sensitivity (log IC50 concentration) of a matching approach in combination with NMF signature extrac- panel of 930 cell lines (separately for 29 cancer types that had a tion, we again used the set of simulated cell line exomes as above, sufficient number of cell lines available) to a set of 518 drugs from where the true somatic mutation signatures are known because the GDSC database37. In particular, we used Random Forest (RF) the matched-normal was available. The ancestry-matching regression applied to the complete set of genetic or epigenetic significantly improves the cosine similarities towards true NMF markers (listed above) in an individual cell line as features signature spectra, compared with the usual approach of filtering (Fig. 2a). In addition to mutational signatures inferred herein, we population variants (p = 0.021, Wilcoxon test; Supplementary also considered the cell line mutational signatures reported by Fig. 1g) and similarly so for the signature exposures (p = 0.017; two recent studies31,32, obtained using approaches that did not Supplementary Fig. 1g). We conclude that our implementation of account for ancestry and that have moreover fit the data to pre- ancestry-matching benefits NMF mutation signature extraction in existing sets of SBS signatures, rather than extracting signatures unmatched cancer samples; we recognize that future variations on de novo from cell line genomes (see Methods section). this methodology might bring improvements. Firstly, mutational signatures predicted drug sensitivity We jointly inferred trinucleotide (or SBS) signatures together significantly better than all other tested types of alterations: with a set of indel mutational features. Examining the SBS part of the increase in accuracy of RF models over CNAs, DNA the spectrum, this yielded 30 cell line mutational signatures that hypermethylation, and oncogenic mutation features is significant very closely match (at a cosine similarity cutoff ≥0.95) the known (at p < 0.05, by corrected Friedman test followed by the post-hoc tumor SBS signatures, and a further 22 cell line signatures that Nemenyi test on ranks; Fig. 2b; Methods section; average rank for match known SBS tumor signatures (at a stringent cosine mutational signatures was 3.75 while for other types of (epi) similarity ≥0.85 and <0.95; a randomization test estimated that genetic features it was 4.20–4.49, considering all RF models). a ≥0.85 cosine threshold corresponds to a 1.8% FDR in matching Secondly, mutational signatures found by our ‘ancestry matching’ the correct SBS, Supplementary Fig. 4b). approach perform significantly better than other cell line The former group was labeled with the name of the signatures recently reported31,32, comparing across the set of cell corresponding SBS signature, while the latter similarly so plus lines that overlap between the publications (Fig. 2b, c). Moreover, the suffix “L” (for “like”). In some cases, our cell line signatures these previous sets of cell line mutational signature exposures were similar to more than one previous tumor SBS (Supplemen- were less predictive of drug sensitivity than were CNAs, tary Fig. 4a) and they were named such as to make this evident, oncogenic mutations and DNA methylation, suggesting the for instance our SBS26/12L matches the DNA mismatch repair utility of adjustment for germline spectra contamination prior (MMR) failure signature SBS26 and a possible MMR failure to mutational signature inference (Fig. 2b–d). signature SBS1238 (more similar signature listed first). Note that a Next, we applied a different test that considers the average error comparable degree of ambiguity is also observed among some of in predicting the drug sensitivity profile (relative RMSE of a RF the known tumor SBS mutational signatures (Supplementary model in crossvalidation; Fig. 2a), averaged across all drugs in a Fig. 5). The full set of 52 mutational signatures we inferred given tissue. In most cancer types, the mutational signatures and their ‘exposures’ across cell types are visualized in obtained herein were better predictors of drug response, than all Supplementary Figs. 2 and 3, and corresponding data is provided the usual genomic and epigenomic features (13 out of 16 tested as Supplementary Data S2 and S3. cancer types compared with DNA hypermethylation, 10 out of 16 Additionally, there were five mutational signatures that appeared for CNAs, and 10 out of 16 for oncogenic mutations). Also, in specific to cell lines (SBS-CL), meaning they did not closely match most cancer types (Fig. 2a) our cell line signatures significantly one of the signatures from current tumor catalogs (Supplementary outperformed recent methods to infer mutational signatures31,32,42 Fig. 2 and Supplementary Data S2). These mutational processes naive to germline mutational spectra (in 22 and 25 out of 27 cancer may be evident only in rare tumor types or they may be active types for the two previous methods that used the same set of cell predominantly in cultured cells rather than in tumors. Some might lines, Fig. 2a, and in 11 out of 16 cancer types on a set of originate from the incomplete separation of other signatures overlapping cell lines for a third method (Supplementary Fig. 7a); 4 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE a ● BONE,GBM,LGG,NB,EWING ● MM,CLL,DLBC,LYMP,ALL,LAML b CD All cancer types Connected=not significant Best 2 3 4 5 Worst ● COREAD,STAD,BRCA,OV,UCEC ● SKCM,KIRC,LIHC,MESO,SARC,THCA Not connected = ● HNSC,PAAD,BLCA,ESCA,LUAD,LUSC,SCLC Exp significant difference Ref. 32 This work Muts HypMet Ref. 31 p−value=0.2019 p−value=0.0467 p−value=0.0222 CNAs 1.10 ● ●● 1.10 ● 1.10 ● ● ●● ● c LUAD RRMSE (HypMet) ● ● CD RRMSE (CNAs) ● RRMSE (Muts) ●● ● ● ● ● ● ● ● 1 2 3 4 5 6 1.05 ●● ● ●● ● ● ● 1.05 ● ● ● ● ● 1.05 ●●● ● ● Exp HypMet ● ● 1.00 ● ● This work CNAs Ref. 32 Muts 1.00 Ref. 31 1.00 0.95 ● ● ● CD COREAD 0.95 1.00 1.05 1.10 1.00 1.05 1.10 1.00 1.05 1.10 1 2 3 4 5 RRMSE (MSigs, this work) RRMSE (MSigs, this work) RRMSE (MSigs, this work) Exp HypMet Predict each This work Ref. 32 p−value<0.0001 p−value<0.0001 Muts Ref. 31 drug RRMSE (MSigs, Ref. 32) RRMSE (MSigs, Ref. 31) ● ● CNAs 1.10 ● ● ● Drug Drug 1.10 ● Sig1 Sig2 … … BRCA ● 1 2 CD ●● ●● ●●● Cell1 1.2 0.6 -1.2 2.9 ● ● ● ● ● ●●● ● ● ● Expo ln 1 2 3 4 5 ● ●● ● ● ● ●●● Cell2 6.5 12.3 3.4 4.1 sure IC50 ● ●● ● Cell3 2.5 13.1 -2.2 -4.8 ● ● 1.05 ● ● 1.05 ● … … … ● ● ● Exp Muts ● ● ● Repeat with Evaluate This work Ref. 31 2nd marker set (5X10-fold XV) CNAs Ref. 32 ● HypMet RRMSE_MSig MSE1 MSE2 RR RR … 1.00 RRMSE_CNA avg ... 1.00 1.05 1.10 1.05 1.10 RRMSE (MSigs, this work) RRMSE (MSigs, this work) d Predictive models (%) 50 e exp_only either combination_only other_only Exp Exp+MSigs 25 Exp+CNAs Exp+HypMET 20 Exp+Muts M As yp s ET s gs s+ xp Si p s +C T C s 0 20 40 60 80 H NA ut M Hy ut g gs ME Si ig E Si N M M + M Predictive models (%) Other Combination M MS Fig. 2 Prediction of drug response with mutational signatures and other molecular data types. a Predictive performance (RRMSE, relative root-mean- square error) of drug response prediction with mutational signatures (“MSigs") reported here and previously31,32 and other data types (oncogenic mutations (“Mut”), copy number alterations (“CNAs”) and DNA hypermethylation (“HypMet”)). P-values of paired one-sided Wilcoxon signed rank test are reported on the plots. Dashed line denotes the diagonal. Bottom right panel shows a schematic of how RRMSE for each tissue was estimated, where “Sig” is mutational signature or other marker (CNA etc.), “Cell” is cell line, “XV” is crossvalidation, and “Predict” implies a Random Forest model. b Average rank diagrams for the performance (RRMSE) of the drug response prediction from various sets of markers, using Random Forest. “This work” refers to the mutational signatures inferred here, while “ref. 31.” and “ref. 32.” refers to prior sets of mutational signatures. Each graph shows: the ranking among the different marker sets (those at the left-hand side are the best performing) and the significant differences between pairs of marker sets (if their ranks are at least critical distance (CD) apart, the difference in predictive performance is statistically significant at p < 0.05, by the Nemenyi post-hoc test, two-sided). The groups of marker sets for which there is no significant difference are connected by black lines. c Tests shown separately for four tissues-of-origin with the highest number of cell lines in our panels. d The percentage of the Random Forest models that are predictive of drug response, defined as having a predictive error lower than the one of an uninformative default model (predicts the average log IC50 for every cell line). Expressed relative to the total number of testable drug-tissue pairs. e The percentage of models that are predictive with gene expression but not another feature type (“exp_only”), by other feature type but not with gene expression (“other_only”), by either model (“either”), or only by a combination of both (“combination”). Source data are provided as a Source Data file. p < 0.0001, p < 0.0001, and p = 0.041, respectively, Wilcoxon test percentage of predictive RF models when markers are combined for decrease of relative RMSE). than with gene expression alone (Fig. 2d, e). If gene expression Gene expression was overall very highly predictive of drug markers are unavailable, the mutational signatures were still response, (Fig. 2b–d), consistent with recent reports30,41. We complementary to oncogenic mutations, CNAs, or DNA asked if the predictive power of gene expression can be methylation (Fig. 2d). complemented by additionally including mutational signatures Next, we considered the complementarity analyses at the level and/or various sets of genetic markers. We predicted the drug of individual drugs, asking if the profile of drug sensitivity that response profile in RF models as above (Fig. 2a), but here by can be predicted by a combined RF model (e.g. gene expression using combinations of marker types with gene expression, and and mutational signatures) could also have been predicted by the tallying the predictive RF models (drug-tissue pairs with better- two RF models drawing on the individual sets of features – on than-baseline RRMSE in crossvalidation; Methods section). gene expression only, or on mutational signatures only (Fig. 2e). Notably, gene expression is complemented by mutational The number of RF models where gene expression by itself is not signatures and also by other types of features, yielding a higher predictive but mutational signatures are predictive is substantial NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 (891 drug-tissue pairs, plus 432 where only a combination of copy-number changes (55), or DNA methylation (64) in the same signatures and expression is predictive). This is higher than the test. We list the associations in Supplementary Data S5. number of drug-tissue pairs where gene expression is not Given that the amount of cell lines available to the replication predictive but driver mutations are (563 plus 352), and similarly analysis is reduced and thus statistical power is limiting, so for CNA (633 plus 347). In addition to gene expression, using particularly for some tissues (Supplementary Fig. 13b), we suggest DNA methylation as a baseline also supports that response that some associations at permissive thresholds (here, nominal profiles for many drug-tissue combinations can be predicted p < 0.005) might be of interest for use as supporting evidence, only by mutational signatures (Supplementary Fig. 7b). This corroborating other associations (see below). suggests that the predictive signal in mutational signatures Secondly, we performed an internal replication analysis within does not simply reflect cancer subtype or cell-of-origin, at least GDSC, enforcing that associations must be detected with two or to the extent that subtype can be identified via gene expression or more drugs that share the same molecular target. In total, 228 DNA methylation patterns. Overall, mutational signatures, drugs in the GDSC data could be tested in this “same-target” considered collectively, can complement gene expression and analysis. Effectively, multiple drugs serve as pseudoreplicates, and other types of markers in predicting drug sensitivity profiles of additionally this test may help discard associations due to off- cancer cells. target effects, which are more likely to differ between two drugs than their on-target effects. Here, we identified 971 significant associations for mutational signatures, 206 for driver mutations, Associations with drug response that replicate in independent 288 for copy-number changes, and 762 for promoter DNA data sets. Because of reproducibility concerns in large-scale drug methylation (at effect size d > 1 and FDR < 15%) (Fig. 4a, b; data screens43–45 that might stem from technical reasons or from in Supplementary Data S7 and S8 for the associations between the cancer cell line evolution during culture, we asked if the asso- default FDR<15% threshold, and the permissive threshold at ciations between mutational signatures and drug responses p < 0.005, respectively). replicate across data sets. To this end, we tested for associations Some associations overlapped between the same-target and the involving various markers, with the additional condition that GDSC-PRISM replication analyses, suggesting more robustness: the associations also replicate in an independent data set. We Supplementary Data S8 contains those replicated associations implemented a randomization-based procedure that tests that the that were seen either across multiple cancer types, and/or across smallest effect size (Cohen’s d statistic) across both datasets is multiple drugs that target the same pathway, and/or across above chance (Fig. 3a; Methods section). A value of d ≥ 1, typi- different replication methods (including also CRISPR genetic cally considered a large effect size, implies that a difference in screens, see below). We suggest that this ‘silver set’ of 3911 mean drug sensitivity (log IC50) between the cell lines positive for associations, where each replication at a FDR < 25% is further a marker and those negative for a marker is greater than the supported by one or more additional replications at a suggestive pooled standard deviation (of the log IC50) of the two sets of cell (nominal p < 0.005) threshold, to be potentially suitable for lines. These replication tests were performed on binarized further analyses or follow-up. mutational signatures i.e. the signature present/absent indicator variables (Supplementary Fig. 15b), considering different cancer types individually (Supplementary Data S4). Additionally the Integrating drug screening data with genetic screening data to same tests were performed on the usual markers in cell line obtain robust associations. As a third type of replication ana- screening analyses, including oncogenic driver mutations lysis, we prioritized cancer vulnerabilities by intersecting the drug (Fig. 3b–e), CNAs (Fig. 3e), and promoter DNA methylation30. sensitivity data with genetic screening data sets. Our rationale was We considered three different types of replication analyses: an that a biological process may be targeted similarly by pharma- external replication in an independent drug screening data set, cological inhibition of a protein, or by editing the genes that internal replication with multiple drugs affecting the same target, encode the corresponding protein. In specific, we examined and an external replication using CRISPR/Cas9 gene knockout sensitivity to CRISPR/Cas9-mediated knockouts that target the fitness screening data. Randomization p-values from these three protein-coding genes, across a panel of 517 cell lines48 that replication methods (across various tissues) were rarely inflated overlapped the GDSC cell lines. We identified many associations for mutational signatures, and in fact commonly exhibited (at effect size d > 0.5 and FDR<15%) with conventional markers – deflation (mean lambda across tissues 0.68–1 for different oncogenic driver mutations (n = 123), copy number alterations replication methods, Supplementary Fig. 8) suggesting an over- (n = 100), and DNA methylation (n = 86) – that replicated across all conservative bias in the replication test as implemented. drug and genetic data (Supplementary Data S6, Figs. 4 and 3e, The few tissue-method combinations that did exhibit inflation in and Supplementary Fig. 9b, c). p-values (lambda >1.3; Supplementary Fig. 8) were omitted from Demonstrating the utility of this test, we recovered well-known further analyses of the associations; the full set of associations examples of the oncogene addiction paradigm, such as the breast are nonetheless included in the Supplementary data, for and esophageal/gastric cell lines with ERBB2 (HER2) amplifica- completeness. tion being sensitive to inhibitors of EGFR and ERBB2 (afatinib Firstly, we performed a replication analysis where the drug and 4 other drugs), and also to the ERBB2 gene knockout association data from the GDSC was tested against another drug (replication FDRs all <15%). Similarly, we recapitulated the screening data set: PRISM (derived from an experimental known associations between amplifications of a chromosomal methodology based on pooled barcoded screens46; 348 cell lines segment 7q31 including the MET oncogene, and sensitivity of and 178 drugs overlap with the GDSC set). In total, 290 drug- esophageal/gastric cancer to crizotinib49,50 and also 3 other MET mutation signature associations were robustly supported across inhibitors (FDRs < 6%; Fig. 3E). The BRAF mutations in skin and both GDSC and PRISM (d ≥ 0.5 and same direction of effect in colorectal cancer likewise sensitize to different BRAF inhibitors in both datasets and additionally requiring randomization test both the GDSC and PRISM drug data, and also to BRAF gene FDR<15%; adjustment using the q-value method47), observed disruption (Fig. 3b and Supplementary Fig. 9d). Conversely, we across diverse tissues and diverse signatures (Fig. 4a, b and note that some oncogene mutations can also confer drug Supplementary Fig. 12a, b). This exceeds the number of drug resistance e.g. NRAS-mutant leukemia cells (Supplementary associations replicated in PRISM involving driver mutations (37), Fig. 10, n = 38 hits in Fig. 4a), consistent with prior reports 6 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE Fig. 3 Detecting robust associations between genetic or epigenetic markers and drug sensitivity by replication across measurements. a A schematic of the randomization test methodology to detect replicating associations using three different tests: (i) consistent effects of a drug between two screening assays (GDSC/PRISM), (ii) effects of a drug consistent with effects of knockout of the target gene (GDSC/PSCORE), and (iii) effects consistent across different drugs that share the same molecular target (GDSC/GDSC or “same target”). Box plots and scatterplot are illustrative. b–d Examples of replicated associations of a known example of oncogene addiction (to BRAF, b) and of additional cancer vulnerabilities associated with mutations in tumor suppressor genes (ARID1A, c; TP53 d). Y-axes show a Z-score derived from either the ln IC50 value (i.e. drug sensitivity, in columns labeled “GDSC” or “PRISM”) or from the CRISPR essentiality score (in the column labeled “PSCORE”). Horizontal brackets show FDR for replicated significant difference between wild-type and mutant genotypes, obtained via a randomization test in panel a, where color denotes the type of the replication test (“GDSC”, “PRISM” or “PSCORE”). The center line of box plots denotes medians and the hinges correspond to the 1st and 3rd quartiles, while whiskers extend to 1.5× IQR from the hinges. e, Effect sizes of markers that associate with drug response in the GDSC drug screen (X-axes) and with response drug target gene knock-out in the Project SCORE genetic screen (Y-axes), shown separately for copy number alterations (CNA) and mutations in cancer genes (Muts). Gray points represent all tested associations, while colored points denote the statistically significant associations that also meet an effect size threshold. Blue lines are the contours of the 2D kernel density estimates. Representative points are labeled. Drug names on grouped labels (“Muts” sub-panel) are ordered by their appearance on the plot from left to right. Source data are provided as a Source Data file. NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 a GDSC/GDSC (same target) GDSC/PRISM GDSC/PSCORE Signatures Muts CNAs SBS43/60L SMAD4 loss:LUSC5 (HRAS,IRF7) SBS−CL3 PIK3R1 loss:LUAD7 SBS55L/60L PBRM1 loss:LIHC7 (CDKN2A) SBS3L MLL2 loss:KIRC24 (FAT1,IRF2) SBS49 MAP3K1 loss:GBM110 SBS40L FAT1 loss:DLBC2 SBS−APOBEC CTNNB1 loss:COREAD53 (ARFGAP3) SBS−CL4 CEP290 gain:LUSC32 (BCL11A,REL,XPO1) SBS44L BRWD1 gain:LUSC30 SBS53L/38L STK11 gain:LUSC27 (EIF4A2,RFC4) SBS6 NOTCH1 gain:ESCA11 (CCND1,CTTN) SBS24L NFE2L2 gain:BRCA47 (EGFR) SBS17b CUX1 loss:STAD42 SBS17aL BRCA1 loss:STAD41 SBS13 BMPR2 loss:GBM112 SBS33L RB1 loss:GBM111 SBS7dL NR4A2 loss:DLBC4 (CDKN2A) SBS31L MAP2K4 gain:LUSC48 (MCM3) SBS−CL1 KEAP1 gain:LUSC26 (SOX2) SBS42L AXIN2 gain:BRCA31 SBS10a/56L ASXL1 loss:LUAD31 SBS29L/18L ARID2 loss:COREAD28 (SOX9) SBS2 SMARCA4 loss:BRCA33 SBS1 CREBBP gain:GBM26 (MDM2,NUP107) SBS16L CTCF gain:GBM122 SBS52 CDH1 gain:ESCA10 SBS46L SOS2 gain:COREAD47 (MYC) SBS28 KDM6A loss:STAD23 SBS26/12L PTEN loss:LUSC8 SBS11 EP300 loss:LUSC16 SBS45/38L/4L ATM gain:BRCA27 (CLTC,PPM1D) SBS36/56L MLL3 gain:BRCA10 (MDM2,NUP107) SBS23/11L BRAF gain:BLCA19 (ERCC5,ING1,IRS2,TFDP1) SBS38/45L MAP4K3 loss:LUSC22 (SMARCA4) SBS48L ARID1B gain:STAD45 (MET) SBS22L ARID1A gain:LUAD32 SBS4/45L PIK3CA loss:KIRC23 SBS20 ELF3 gain:GBM61 (MDM4,PIK3C2B,ZC3H11A) SBS59 BPTF gain:BLCA9 SBS60 KRAS loss:LUSC10 (CDKN2A) SBS54 EGFR gain:BLCA15 (ASXL1) SBS8L/4L NRAS gain:BRCA26 (CDK12,ERBB2,MED24) SBS18/36L MGA loss:LGG7 (RPL22) SBS25L TP53 gain:LGG16 (EGFR) Other Other Other 100 50 0 50 20 0 20 40 100 0 100 b Number of associations (Sensitivity | Resistance) GDSC/GDSC (same target) GDSC/PRISM GDSC/PSCORE n signatures (FDR<15%) 20 15 10 5 0 SCLC BRCA SARC LUAD DLBC MESO PAAD SKCM MM HNSC GBM/LGG LIHC LUSC LAML ALL/CLL COREAD_MSS ESCA/STAD EWING NB Fig. 4 Tally of the significantly replicated associations of drug sensitivity or resistance with mutation signatures and other markers. a Comparison of the number of statistically significant associations (FDR<15% by randomization test, additionally requiring an effect size d > 0.5 for GDSC/PRISM and GDSC/PSCORE tests, and d > 1 for the GDSC/GDSC (same-target) test) per feature, among mutational signatures (“Signatures”), oncogenic mutations (“Muts”) and copy number alterations (“CNAs”) in the three types of replication tests (see Fig. 3a). Features are ranked by the total number of significant associations, either for drug sensitivity (negative side of X-axis) or resistance (positive side of X-axis). b The number of different mutational signatures that have statistically significant associations across various cancer types (at FDR<15%; we consider signatures that have >1 significant association per cancer type), in the three replication tests. Source data are provided as a Source Data file. (discussed in Supplementary Note 1). These and other striking biological roles of these genes, as well as prior preclinical studies associations with gene mutations, CNA, and promoter DNA (details in Supplementary Note 1). This demonstrates the power methylation that replicated in the genetic screening data are of the joint analyses of drug and genetic screening data here and highlighted in the global data overview in Fig. 3e and elsewhere51, suggesting that the other associations we identified Supplementary Fig. 9c. here (Supplementary Data S6) provide cancer dependencies In addition to oncogene addiction, replicated associations promising for follow-up. For example, our analysis identifies can suggest ways to target mutated tumor suppressor genes via vulnerabilities of TP53-mutant cells to manipulating the activity synthetic lethality. One example is a CDKN2A deletion that of the CHK2/CDC25A/CDK2 axis across five different cancer sensitizes brain cancer cells to palbociclib (a CDK4/6 inhibitor) types (Fig. 3d and Supplementary Fig. 9e), echoing prior work on and to knockouts in CDK4 and CDK6 genes. In contrast, RB1 therapeutic interventions on CHKs in TP53-deficient cells mutations were associated with resistance, consistent with the (Supplementary Note 1). These examples also illustrate how an 8 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE integrated analysis of drug screening data with genetic screening for WRN activity. Overall, the ability to recover the known WRN data can reveal drug effects exerted via secondary drug targets dependencies of MMR-deficient cell lines estimated via trinu- (e.g. likely CHK2 for the MK-8776 inhibitor of CHK1; Fig. 3d, see cleotide mutational signatures supports the utility of our meth- discussion in Supplementary Note 1). odology to infer mutational signatures in cell line genomes. We highlight a robustly supported synthetic lethality example Beyond WRN disruption, the MMR signatures as well as other involving mutations in the ARID1A tumor suppressor gene and mutation signatures can predict sensitivity to many perturbations, the inhibition of the AKT2 gene or protein (Fig. 3c). In particular, including those that the MSI status nor the other genetic markers ARID1A mutant colorectal cell lines are more sensitive to the would not predict (Supplementary Fig. 9a; we note that the knock-out of the AKT2 gene by CRISPR, as well as to the pan- converse is also true, at least for the few cancer types where MSI AKT inhibitors GSK690693 and capivasertib/AZD5363 (FDR = labels are available). 6% and 12% in the replication test, respectively). The same is Next, we systematically examined all mutational signatures for observed in ovarian cancer cell lines, again involving AKT2 the overlap between associations in the GDSC drug screen and knockout and the same two inhibitors (at FDR = 9% and 12%, Project SCORE genetic screen. This yielded 130 associations (at a respectively). This is supported by additional AKT inhibitor randomization FDR < 15%, and additionally requiring an effect drugs: afuresertib (FDR=6%), AKT inhibitor VIII (FDR = 21%), size of d > 0.5 in both the genetic and the drug screens) that and uprosertib (FDR = 5%) in colon, and MK-2206 (FDR = 9%) replicated across data sets – a higher number than for oncogenic in ovary (Supplementary Data S6). Further evidence for an driver mutations, CNAs, and DNA methylation (123, 100, and 86, interaction between these genes is found in tumor genomic respectively, at the same FDR threshold). These associations analysis. The AKT2 oncogene can be amplified in ovarian, (Fig. 4, Supplementary Fig. 12c; full list in Supplementary endometrial, pancreatic and other cancer types, while the ARID1 Data S6) involved k.o. in 64 different genes, indicating that tumor suppressor commonly bears truncating mutations in many mutational signatures associate with a variety of target genes and cancers. In tumor genomes, AKT2 alterations significantly co- suggesting potential points of intervention for follow-up. The occur with ARID1A alterations (OR = 2.0, FDR<0.1% in MSK- number of replicated associations involving mutational signatures IMPACT cohort of 10,945 samples;52 replicated at OR = 1.4, highly ranked by this replication analysis (such as the FDR<0.1% in an independent TCGA pan-cancer cohort of chemotherapy-associated SBS25L, the haloalkane exposure- 10,967 samples; analysis via cBioPortal53). These genomic associated SBS42L, and the signature possibly related to NER associations support that the AKT2 amplifications may bring a deficiency61 SBS8L/4L; Fig. 4b, Supplementary Fig. 12a; n = 12, 8, selective benefit to ARID1A-mutant tumors. Overall, our analyses and 7 replicated associations at FDR 15%, respectively) broadly solidify the notion that the PI3K/AKT/MTOR signaling inhibi- matches the number of replicated associations involving common tion is a vulnerability of ARID1A-mutant cells54–57, as reported driver mutations such as EGFR or TP53 or KRAS (n = 18, 17, and before for individual examples of cell lines sensitized to AKTi 13, respectively), or known copy number change events such as drugs upon silencing of ARID1A56, and we further suggest ERBB2 gain (n = 17 replicated associations in Project SCORE) specifically AKT2 as an opportune point of intervention. (Fig. 4a and Supplementary Data S6). We also show the tally of Next, we applied this same statistical methodology (Fig. 3a; associations at a more permissive 25% FDR in Supplementary Methods section) to identify replicated drug sensitivity associa- Fig. 12, further supporting how mutational signatures provide tions involving mutational signatures in cell line genomes. markers as commonly associated with drug response as the usual markers based on driver mutations and CNA. We note that in this and further analyses we have, Mutational signatures associated with sensitivity both to conservatively, stratified colorectal cancer cell lines into MSI pharmacological and genetic perturbations. As a positive con- and MSS (association counts without stratification are in trol in a study of mutational signatures as markers, we considered Supplementary Fig. 13a and Supplementary Data S9), based on a recently reported vulnerability of cell lines that are MSI being common in colorectal cell lines, and MSI status being microsatellite-instable (MSI) and therefore deficient in the DNA strongly associated with mutational signatures9,38,62. mismatch repair (MMR) pathway, which do not tolerate the loss of the WRN gene48,58–60. MMR deficiencies in tumors are known to associate with MSI and with trinucleotide mutational sig- Robust associations involving mutational signatures replicate natures SBS6, 15, 21, 26, and 449 (and additionally SBS14 and 20 across multiple cancer types. We next focused on those drug which result from MMR failure concurrent with deficiencies in associations involving mutational signatures that recurrently replicative DNA polymerases). In a joint analysis of MSI-prone replicated across more than one replication method (see Fig. 3a) cancer types (colorectal, ovary, stomach, uterus), we found links and/or more than one cancer type (‘silver set’, Supplementary between the MMR SBS signatures that we inferred in the cell line Data S10). We noted that some associations in this set recurred in exomes and the sensitivity to WRN knockout. However, levels of three or more methods and/or tissues, thus we introduced a more statistical support were variable across the MMR signatures stringent tier of hits, with a higher priority for follow-up. These (FDRs <0.01%,<0.01%, <0.01%, 11%, 18%, 21%, and n.s. ‘golden set’ hits recurred in at least three different tissues or in [associated with resistance]) for SBS20, 15, 6, 26/12L, 14, 21, and three different tests with effect size d > 0.5 and significant at 44L, respectively; Supplementary Fig. 11a). Additionally, we p < 0.005, and at least once at FDR < 25%. This resulted in 995 noted some additional signatures with a high weight on the indel higher-priority associations (tallying both mutational signature components – and thus might be MMR-related: SBS33L, SBS54, associations, and the driver mutation and CNA associations; and SBS-CL1 (Supplementary Fig. 2) – also predicted sensitivity Supplementary Data S8). to WRN loss (Supplementary Fig. 11), in case of SBS33L with a A common occurrence in this higher-confidence association high effect size. Thus, some MMR-associated signatures are more set was involvement of mutational signatures that were associated robust markers for WRN inhibition (particularly the C > T-rich with DNA repair failures in previous analyses of tumor genomes. SBS15 and SBS6, as well as SBS20) than the other MMR failure- This included: MMR failures (various SBS; listed in Fig. 5a, b), associated signatures (such as the T > C rich 26, 21, or 44). BER failures (SBS36/56L15,63,64, SBS30L/7bL/11L14,65), likely Conceivably, this might be because these signatures reflect dif- NER failures (SBS8L/4L) and replicative DNA polymerase failures ferent types of MMR failure that confer differential requirements (particularly SBS14 and SBS20; additionally SBS56/10aL/36L may NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 a 5 AKTi SBS14 SBS20 SBS21 SBS6 b AKTi Sensitivity Resistance Number of associations associations (Resistance | Sensitivity) associations SBS15 SBS26/12L SBS44L (pvalue < 0.005) 2.5 40 −log(pvalue) 30 0 20 10 0 2.5 L 5 S6 4 4L 0 12 S1 S1 S2 S4 SB 6/ SB SB SB SB S2 5 SB c * * Sensitivity Resistance * DNA Signature previously associated with chemical exposure Number of associations 150 ** * ° repair signatures * (FDR<25%) 100 °° * * ° * ° ° 50 ° ° °°° °° 0 SBS56/10aL/36L SBS30L/7bL/11L SBS−APOBEC SBS45/38L/4L SBS53L/38L SBS29L/18L SBS10a/56L SBS55L/60L SBS−MMR1 SBS−MMR2 SBS43/60L SBS26/12L SBS12/26L SBS18/36L SBS38/45L SBS36/56L SBS23/11L SBS4/45L SBS8L/4L SBS−CL1 SBS−CL3 SBS−CL4 SBS−CL2 SBS−CL5 SBS17aL SBS17b SBS44L SBS33L SBS40L SBS19L SBS10b SBS25L SBS42L SBS22L SBS24L SBS16L SBS48L SBS31L SBS7dL SBS32L SBS46L SBS14 SBS3L SBS13 SBS49 SBS7a SBS5L SBS7b SBS21 SBS11 SBS60 SBS54 SBS59 SBS52 SBS20 SBS28 SBS15 SBS2 SBS6 SBS1 d associations(pvalue < 0.005) Enrichment of resistance 0.6 log2 fold change 0.4 0.2 0.0 ERK MAPK signaling Chromatin histone Chromatin histone and degradation EGFR signaling DNA replication Protein stability Other, kinases ABL signaling JNK and p38 p53 pathway Unclassified PI3K/MTOR methylation acetylation signaling signaling Fig. 5 Highlighted examples of robustly supported associations involving mutational signatures. a All tested associations between AKT inhibitors and DNA mismatch repair mutational signatures, across all three replication tests (by the “two-way” randomization test, see Methods section). Each bar represents a p-value of one association. For associations with effect size <0.2, p-values were not calculated in the randomization procedure and are here shown as having p = 0.5. b Associations having p < 0.005 between AKT inhibitors and individual DNA mismatch repair signatures. c The tally of significant associations at FDR < 25% across all three replication tests. The * and ° symbols denote groups of mutational signatures; see key embedded within the panel. d The average log2 ratio of observed vs. expected frequencies of occurrence of drug target pathways in resistance associations (p-value < 0.005) with signatures of chemical exposure, over the top 6 signatures by the number of resistance associations (SBS25L, SBS18/36L, SBS42L, SBS11, SBS22L, and SBS4/45L). Source data are provided as a Source Data file. be in this group). As tentatively DNA repair-associated By analogy to this, we also identified two opposing trends in signatures, here we additionally considered6,9,66 SBS18/36L based the mutational signatures association tally. Firstly, the signatures on the similarity of the spectrum to the SBS36/18 and because it associated with DNA repair failures, as listed above, tend to be was found in MUTYH-variant patients63,67, and additionally more often sensitizing (considering relative frequencies of SBS33L, SBS54, and SBS-CL1 because they have prominent indel sensitivity to resistance associations shown in Fig. 5c; also see components and were associated with sensitivity to WRN loss top of the plot in Supplementary Fig. 12a for breakdown by type (Supplementary Fig. 11). Those DNA repair-associated signatures of test). encompass 278 of 701 associations involving mutational signa- Secondly, there is a group of mutational signatures tending tures in this high-priority set; some individual examples are towards resistance associations; these signatures also exhibit a discussed below. Therefore, mutational signatures resulting from higher overall number of associations (Fig. 5c; bottom of plot in DNA repair failures often result in drug vulnerabilities. Supplementary Fig. 12a). In this group, 6 out of top-7 signatures When inspecting the overall balance of sensitivity versus were previously linked to exposures of mutagenic chemicals: resistance associations, we noted that driver mutations and CNA SBS25L (unspecified chemotherapy), SBS18/36L (reactive oxygen present a mix of sensitizing and resistance associations. These species), SBS42L (haloalkane exposure68), SBS11 (a DNA may be unevenly distributed across genes: see NRAS example methylating agent, the drug TMZ), SBS22L (an agent generating mentioned above, biased towards resistance, while EGFR is biased bulky DNA adducts, aristolochic acid69), and SBS4/45L (mix of towards sensitivity (perhaps because mutant EGFR is a target for agents from tobacco smoke, where the mutagenesis likely results many approved drugs). mainly from DNA adducts by polycyclic aromatic hydrocarbons). 10 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE Of note, among these, the oxidative damage SBS18/36L does Mutational signatures associated with various DNA repair show some sensitizing associations as well (Fig. 5c). These six failures predict drug sensitivity. A manual curation of the sets of signatures of mutagen exposures associated with resistance to robust associations (Supplementary Data S8) reveals that several various drugs, which are overall enriched in drugs targeting e.g. MMR signatures associate with sensitivity to AKT serine/threo- chromatin histone modification, DNA replication, the p53 nine kinase inhibitors (Fig. 5a, b). This is seen consistently across pathway, JNK and p38 signaling, and ABL signaling (Fig. 5d many tissues: in colorectal, skin, lung (small cell), and brain and Supplementary Fig. 15 breakdown per signature). This data (associations at FDR<15%), and additionally in prostate, ovary, suggests that, overall, mutational signatures of prior chemical and stomach/esophagus cancers (associations at permissive FDR exposures in cancer cells commonly predict resistance to future thresholds, all with p < 0.005); see Fig. 6a–c for examples invol- drug exposures. ving SBS26/12L, SBS14, and SBS20, respectively. The associations GDSC/GDSC (same target) GDSC/PRISM GDSC/PSCORE a SBS26/12L ● Absent ● Present ln IC50 or CRISPR score (standardized) COREAD COREAD_MSI ESCA/STAD OV 5.0 0% 14% 5% 31% 5% 33% 17% 11% 51% 53% 31% ● ● 49% ● 53% 31% 2.5 ● ● 51% ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ●● ● ●● ● ●● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ●●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ●● ●● ● ●● 0.0 ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●●●●●●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ●● ● ●● ● ●● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ●● ●● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● −2.5 ● ● ● ● ● ● (PSCORE) (PSCORE) (PSCORE) (PSCORE) Afuresertib Afuresertib Afuresertib Ipatasertib Ipatasertib Uprosertib Uprosertib MK−2206 MK−2206 AZD5363 AZD5363 AZD5363 (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) AT7867 AKT2 AKT1 AKT2 AKT3 SBS−CL1 ● Absent ● Present b SBS14 ● Absent ● Present c SBS20 ● Absent ● Present d LAML LUAD NB ln IC50 or CRISPR score (standardized) ln IC50 or CRISPR score (standardized) COREAD_MSS PRAD SKCM COREAD GBM/LGG SCLC 6% 4 19% 4 5% 56% 5% 12% ● 2% 5% 11% 37% 19% 12% 53% 2 ● ● ● ● 3 34% ● ● 20% ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 34% 37% ● 2 ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ●●●● ● ● ● 2 ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ●● ●● ● ●● ●●●● ● ● ●● ● ● ● ● ● ● ●● ● ●● ●● ●●● ●●●●● ●● ● ● ●● ● ● ●● ●● ● ● ● ● ● ● ●● ● ● ●●● ●● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ●● ●● ● ●● ● ●● ●●●● ●● ●●● ● ●● 0 ●● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ●●● ● ●● ● ● ● 1 ● ● ● ●●● ● ●● ● ● ● ● ● ● ● 0 ● ● ● ●● ● ●● ● ●● ● ●● ●● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ●●● ●●● ● ● ●● ●●●● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ●● ● ● ●● ● ●● ● ● ● ●●●●●● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●● ● ● ● ●●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ●●● ●● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● 0 ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ● ● −2 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● −2 ● ● ● ● ● ● ● ● ● ● −1 ● ● ● ● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● −4 ● ● ● ● GSK690693 GSK690693 GSK690693 AKT inh.VIII (PSCORE) (PSCORE) A−443654 A−443654 (PSCORE) (PSCORE) Afuresertib Afuresertib Rucaparib Rucaparib Afuresertib Ipatasertib Uprosertib Uprosertib Uprosertib Rucaparib Veliparib (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) Olaparib Olaparib (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) PARP1 AKT2 AKT2 AKT2 SBS18/36L ● Absent ● Present SBS−CL1 ● Absent ● Present e f ALL/CLL LUSC PAAD SARC BRCA COREAD MESO PAAD 5.0 30% 8% 31% 4 17% ● 20% 7% 11% 35% 7% ● 53% ● 5% ● 46% 46% 2.5 28% ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● 2 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●● ●● ●● ● ● ●●● ●● ● ● ●● ●● ● ●● ●● ● ● ●●● ●● ● ● ●● ●● ● ●● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ●● ● ● ●● ● ●● ● ● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●● 0.0 ●● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ●● ●●● ● ● ● ●●● ● ●● ● ●● ● ●●●● ● ● ●●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●● ●● ●● ● ●● ●● ● ●● ●● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●●● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● 0 ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ●●●● ● ● ● ● ● ● ●●● ● ●● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ●●● ●● ● ● ● ●● ●● ● ● ● ● ● ●● ●● ● ● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2.5 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 ● ● ● ● ● ● (PSCORE) (PSCORE) (PSCORE) (PSCORE) AZD6738 (PSCORE) Tenovin−6 (PSCORE) Tenovin−6 (PRISM) Selisistat VE−822 VE−822 AZD6738 AZD6738 AZD6738 (GDSC) (GDSC) Selisistat Selisistat Selisistat (GDSC) (GDSC) VE−822 (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) (GDSC) VE821 VE821 SIRT1 SIRT1 SIRT1 ATR ATR ATR Fig. 6 Associations with drug sensitivity or resistance that replicate in independent datasets. a–f Examples of associations of mutational signatures with drug sensitivity that replicated (using tests in Fig. 3a) multiple times, across different cancer types and/or different types of replication tests. Y-axes show a Z-score derived from either the ln IC50 value (drug sensitivity: “GDSC” or “PRISM” columns) or from the CRISPR essentiality score (“PSCORE” columns). Horizontal brackets show FDR for replicated associations with the presence/absence of a given mutational signature, obtained via a randomization test (Fig. 3a), where color denotes the type of the test (see legend at top of plot). The center lines of box plots denote medians and the box hinges correspond to the 1st and 3rd quartiles, while whiskers extend to 1.5 × IQR from the hinges. Source data are provided as a Source Data file. NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 involve 9 different AKTi drugs including uprosertib, MK-2206, permissive FDR thresholds (with notable effect sizes, however; and ipatasertib. Several of these drugs have undergone clinical Supplementary Fig. 14a, b). trials showing varying outcomes in unselected patients70,71, Some associations were found involving mutational signatures highlighting the need for identifying predictive biomarkers of recovered from cell line genomes that did not closely match a response to AKT inhibitors. We considered the possibility that known SBS spectrum from tumor genomes. An interesting different MMR signatures have varied utility as AKTi markers; example are associations of the SBS-CL1, an indel-rich signature indeed, the MMR signature SBS26/12L most commonly asso- seen in various cancer types (Supplementary Figs. 2 and 3). Our ciated with sensitivity across different AKTi drugs, with a lower analysis suggests this associates with sensitivity to DNA damage utility of other signatures (Fig. 5a, b). These associations between signaling drugs. Firstly, we identified associations of exposure of MMR signatures and AKTi sensitivity may be mechanistically SBS-CL1 with sensitivity toward PARP inhibitors olaparib, related to associations between ARID1A mutations and AKTi rucaparib, veliparib, and PARP1 gene k.o., observed across three sensitivity that we described above (Fig. 3c). Such a link would be cancer types (Fig. 6d). Secondly, we identified associations with consistent with a reported loss of MMR activity in ARID1A- sensitivity to ATR inhibitors AZD6738, VE-822, VE-821 or k.o. mutant cells72, and with correlations between ARID1A loss in of the ATR gene in four cancer types (Fig. 6f). This example tumors and MMR deficiencies reported in multiple cancer suggests the utility of indel signatures in predicting drug response, types73–75. in this case to a category of DNA repair drugs that are trialed in In addition to MMR, signatures resulting from failures in other the clinic, the ATRi. DNA repair pathways may yield sensitivity associations (Fig. 6d–f Additional associations were recurrently observed across and Supplementary Fig. 14). An example is the signature SBS36/ multiple cancer types. Some examples highlighted by a manual 56L, possibly indicating failed BER since SBS36 was previously curation of these ‘golden set’ associations include: the SBS17aL associated with loss-of-function in MUTYH. In five cancer types, signature and the ZSTK474 drug and PIK3CA/PIK3CB genes; SBS36 was associated with sensitivity to inhibition of EGFR or of SBS3L signature and fedratinib and JAK2 gene; SBS8L and ERBB2 via e.g. afatinib or AST-1306 drugs. These associations midostaurin; SBS17b and fludarabine (and possibly more were additionally supported in CRISPR k.o. of the EGFR or generally DNA antimetabolites). These associations are some ERBB2 genes in skin, liver, and head-and-neck cancer (Supple- representatives among many other examples with a similar degree mentary Fig. 14c). The related signature SBS18/36L also predicts of confidence (based on the FDRs, and on recurrence across sensitivity to these agents in three of the five cancer types independent tissues/replication tests) in the ‘golden set’ of (Supplementary Fig. 14d). Overall statistical support for sensitiv- associations (Supplementary Data S8). In addition to our manual ity associations with SBS36/56L and SBS18/36L was higher for curation, we also provide a prioritization based on a pooled the EGFR-targeting drugs than for all other classes of drugs in p-value across tissues and different replication tests to highlight these cancer types (Supplementary Fig. 14e). Prior studies the top 10 associations with mutational signatures, and suggested that EGFR activity can control various DNA repair additionally with driver mutations/CNA/DNA hypermethylation, mechanisms76–78. in Supplementary Data S11. We further highlight an example involving a signature SBS18/ 36L, associated with DNA damage by reactive oxygen species, and possibly also with certain deficiencies in the BER pathway due Discussion to the similarity of the spectrum with signature 36. This signature A classical way to treat tumors is to employ DNA damaging drugs was associated with sensitivity to two inhibitors of sirtuin (SIRT) and ionizing radiation to target the lessened or overwhelmed proteins, selisistat and tenovin-6, in pancreatic adenocarcinoma, capacity for DNA repair in cancer. Since this may manifest as a lung squamous cell carcinoma, sarcoma, and lymphoid leukemia mutator phenotype, we asked if mutational signatures observed in (all tissues at FDRs ≤ 30%; Fig. 6e). In three out of four tissues cancer cells can serve as markers for treatment by drugs or the associations with SIRTi replicated in the CRISPR k.o. by gene editing. This systematic study generalizes over the phenotype of the SIRT1 gene (Fig. 6e). This adds confidence known individual examples of mutational patterns stemming that SIRT1 may be a promising vulnerability of tumor cells that from deficient HR24,26,27 or MMR, which can guide therapeutic are undergoing and/or have previously undergone oxidative strategies58–60. damage to their DNA and/or have lowered ability to repair such Cancer cell line panels that were screened for drug damage (the current analysis does not distinguish between these response37,79 and for gene loss effects48,80 provided a resource to scenarios). test our hypothesis. However, the lack of matched healthy tissues Another example of a sensitizing mutational signature was means that extracting somatic mutational signatures using interesting due to its occurrence across multiple tissues. The existing methods31,32,42 is a challenge. Since germline variation is SBS30L/7bL/11L, which is ambiguous but possibly linked to base abundant compared to somatic mutations, even slight variations excision repair failures (since the SBS30 was previously associated of germline spectra between populations34,36 can affect the tri- with NTHL1 loss-of-function14,65), associates with sensitivity to nucleotide context mutation tally. We thus subtracted the two related classes of drugs (Supplementary Fig. 14a, b) that expected germline spectrum given the ancestry, and further converge onto the cytoskeleton. Firstly, there are inhibitors of integrated indel features into mutational signature inference from Aurora kinase A, a protein regulating mitotic spindle assembly the cell line WES. Future refinements in the methodology, as well and stability, including the drugs ZM447439, Genentech Cpd10, as availability of WGS of cancer cell lines, will improve its and GSK1070916 (an Aurora B/C inhibitor with some Aurora A accuracy. For instance, this will permit a more detailed set of activity). These associations are also replicated in the k.o. of the indel descriptors, as applied in the recent tumor WGS mutation AURKA gene (Supplementary Fig. 14a). Secondly, this signature signature studies9, in contrast to the limited set of only four indel SBS30L/7bL/11L associates with sensitivity to the vinca alkaloids features we were able to apply in our WES study. vinblastine, vincristine, and vinorelbine that interfere with Among mutational signatures that we inferred, the number of assembly of microtubules and forming of the mitotic spindle associations with drug response was comparable to that of other (Supplementary Fig. 14b). These associations were observed types of commonly used genomic markers. Hundreds of such across AML and CML leukemia and liver cancer at FDR ≤ 30%, as associations significantly replicated in independent data and well as in colorectal cancer and multiple myeloma at more across multiple tissues. Thus mutational signatures appear to be 12 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE similarly robust predictors as driver mutations or CNA markers. an association analysis such as ours. We note that a related issue We note that the associations we identified could not have could affect also the more established markers i.e. driver muta- resulted from tissue-specific variation in drug response, since tions or CNAs in cancer genes: given the rapid accumulation of tissues were considered individually in the association analyses genetic changes in cultured cancer cell lines82,83, and prevalent (the only exception was merging of some cancer types known to epistasis in cancer84,85, it is plausible that recently occurring, be similar genomically, such as esophagus with stomach cancer, unobserved mutations or CNAs affect the ability to identify drug and glioma with glioblastoma). sensitivity markers from analyses of cell line screening data. An important caveat to our study is that the discovered asso- An interesting observation about mutational signatures asso- ciations do not necessarily imply a causal relationship: the ciated with drug activity is that some likely result not from DNA mutational signature-generating process and the drug phenotype repair deficiencies, but instead from exposure to mutagenic may be only indirectly associated. In other words, the mutational agents. Some of these signatures presented many associations in signature might serve as a marker for another alteration (e.g. a our analyses and were, overall, more commonly associated with driver mutation), which may be the proximal cause of the drug drug resistance rather than sensitivity (Supplementary Fig. 12 and sensitivity or resistance. This is similarly possible with mutational Figs. 4 and 6a, b). For example, this includes SBS25 and SBS11, signatures as with other genetic markers (mutations or CNA in reported to be associated with chemotherapy, or signatures linked cancer genes), and also with gene expression markers. Because with exposure to chemicals causing DNA adducts (tobacco these various sets of genomic/transcriptomic features can corre- smoking, aristolochic acid), and additionally an SBS17a-like sig- late across tumors, prioritizing the likely causal relations for nature as well (where SBS17 was associated with gastric acid further follow-up is a challenge that remains to be addressed; exposure, possibly via oxidative damage to the nucleotide pool). larger cell line panels will be helpful, as well as use of isogenic cell Even though cell lines are not exposed to these chemical agents panels. during culture and thus the signatures are presumably not ‘active’, Furthermore, another issue are false negatives in identifying drug associations are identified with such signatures (Fig. 4). associations – again similarly so with mutational signatures as One possible explanation may be that mutational signatures of with other markers – due to sparse data, where often only few cell different processes are sometimes sufficiently similar such that lines from a cancer type bear each marker. Thus, statistical power they are not easily ‘unmixed’ only from trinucleotide spectra. For may be limited for many markers in current cell line screens. instance, SBS17 (mostly A>C/T>G changes) might result from Additionally, because drug sensitivity measurements can be noisy varied mechanisms that converge onto the same spectrum, some (as evidenced in less-than-ideal agreement between different of which are from chemotherapy exposure, while others may be screening data sets43–45) replication analyses across diverse endogenous86. This highlights an example where current statis- datasets will be conservatively biased. It bears mentioning that tical methods may not reliably deconvolve underlying biological absence of a significant association in our lists does not imply mechanisms. This might be addressed by the use of additional there is not an association of that marker with that drug, but mutational features, such as penta- or hepta-nucleotide might rather mean that the analysis at current sample sizes may contexts87, small indels9, copy-number changes22 and strand- be underpowered to detect the association (see e.g. Supplemen- specific, or regional mutation rates10,88. tary Fig. 13b). Another explanation may be that prior exposure to a carci- A key question to be addressed in future work is the clinical nogen would select tumor cells with an altered DNA replication/ relevance of the mutational signature drug markers, which we repair state, which continues after the carcinogen is withdrawn, here identified in cancer cell line panel data. Performing asso- thus generating resistance in cancer cells. Indeed, prolonged ciation studies on tumor genomic datasets (for which the clinical exposure to mutagens that are also cytotoxic – as is the case for data about treatments and patient response are sometimes many cancer therapeutics – is likely to select resistant cells, in available) is complicated by the diversity of therapeutic regimes: some cases via altered DNA replication or repair mechanisms. most patients are treated with multiple overlapping sets of drugs For instance, therapy of tumors with temozolomide (associated and possibly radiotherapy, which makes it challenging to identify with SBS11) is known to select for cells that are resistant due to a effects of individual drugs by retrospective analysis. Additionally, MMR-deficiency via loss-of-function of MSH689. It is conceivable the drug assignment may be confounded by demographics and by that also the various epigenetic changes resulting from carcinogen cancer stage/grade or subtype, further complicating analysis. exposure might confer similar properties. In other words, even a Large controlled randomized trials with treatment and control temporary exposure to a mutagen may prime tumor cells for arms, for which the tumor genomic data is also available, would resisting later drug treatment. facilitate identifying various types of genomic markers (muta- tional signatures or otherwise) relevant to drug response and patient survival. Methods Human cancer cell lines and primary tumor data. We obtained WES bam files With respect to mechanistic insight, future improvements in (human reference genome version hg19) of human cancer cell lines (n = 1072; methodology will refine the mutation signature markers and cancer cell lines from Genomics in Drug Sensitivity in Cancer (GDSC)) from clarify the underlying mechanisms. For instance, a known lim- European Genome-phenome Archive (EGA) (ID number: EGAD00001001039) itation of various mutational signature extraction methods – and WES bam files (human reference genome version hg38) of primary tumors (n = 6154) and their matched normal samples (n = 6154) from the TCGA repo- including ours – is the difficulty of discerning the ‘featureless’ sitory at NCI Genomic Data Commons (access via dbGaP accession phs000178). signatures such as SBS3, SBS5, and SBS881. A further issue that We downloaded samples from the following TCGA cohorts: BLCA, BRCA, COAD, merits attention is timing: a genomic analysis of cell lines31 GBM, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, OV, PAAD, READ, suggested that the activity of some mutational processes is vari- STAD, THCA, UCEC. We aligned the human cancer cell line bam files to the human reference genome (version hg38) using the bwa90 software, and sorted and able in time. While a genome sequence reflects a record of indexed them using the samtools software91. We used the Strelka2 software (ver- mutagenic activity in the past, it may or may not reflect current sion 2.8.4)92 to call single nucleotide variants (SNVs) and small insertions and mutagenic activity, which is presumably more relevant for drug deletions. We called SNVs and indels for cell lines, primary tumors and normal sensitivity phenotypes. Occurrence of recently active signatures is samples. In samples where Strelka2 was unable to run, a re-alignment was per- formed using Picard tools (version 2.18.7)93 to convert the bams to FASTQ and, difficult to identify from bulk DNA sequencing from cell culture, following that, the alignment was performed by executing bwa sampe (version as the recent mutations may not rise to sufficient allele fre- 0.7.16a) with default parameters. The resulting bam files were sorted and indexed quencies to be detected, and may conservatively bias the results of using Picard tools. We used SNVs and indels marked as “PASS” in the Strelka2 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 output. We annotated SNVs and indels with minor allele frequencies (MAF) mutation spectrum. The mutation spectrum contains 96 components, which are obtained from the gnomAD database94 (for SNVs and indels that could be found in the counts of six possible mutation types (C>T, C>A, C>G, A>T, A>G, and A>C, gnomAD). considered DNA strand-symmetrically) in 16 possible 5′ and 3′ neighboring nucleotide contexts for each mutation type. For each cell line, from the cell line’s 96 trinucleotide spectrum, we subtracted the median 96 trinucleotide frequencies of its Data for human cancer cell lines. We downloaded drug response data for 518 TCGA ‘ancestry-matched samples’ (i.e. the TCGA normal samples belonging to the drugs from GDSC (Release 8.3; June 2019)37. We used the natural logarithm of the same population cluster as the cell line in the PC analysis of common germline 50% growth inhibition values (IC50) as a measure of the activity of a compound variants, see above). In case the subtraction resulted in a negative count for some of against a given cell line. If for the same drug activities were available from both the contexts, we set it to zero. GDSC1 and GDSC2 versions, we used GDSC1. We downloaded the information about drugs (the list of drugs’ putative targets and target pathways) from the GDSC website (https://www.cancerrxgene.org/). We manually curated the list to correct Insertion and deletion types. In addition to the standard 96 trinucleotide spectra, inconsistencies (Supplementary Data S12). to extract mutational signatures from cell lines (see below) we considered addi- We obtained the following genetic features of cancer cell lines from the GDSC tional features based on small insertions and deletions. We filtered the regions of repository:30 cancer driver genes (Muts), regions of recurrent focal copy number the exome using the steps (i-iv) as described above. Similarly as for SNVs, we alterations (CNAs), hypermethylated informative CpG islands (HypMet), and discarded indels found at MAF>0.001% in the gnomAD population database and microsatellite instability status (MSI). We downloaded ANOVA_input.txt files for filtered the ones that appeared in >5% of cell line samples. The remaining indels 18 cancer types and pan-cancer analysis from https://www.cancerrxgene.org/ were classified into 4, 5, 8, 14, or 32 different indel types considering: the length of gdsc1000/GDSC1000_WebResources/Pharmacogenomic_interactions.html. We insertion or deletion (considering lengths 1, 2, 3–4, 5+), microhomology at dele- downloaded cancer cell line gene expression data (“sanger1018_brainarray tion sites (considering the lengths of microhomology 1 or 2+), and microsatellites _ensemblgene_rma.txt”) from https://www.cancerrxgene.org/. We selected at indel loci (considering repeat sizes of 1 and 2–5+). We benchmarked the cell line expression data for 1564 genes corresponding to the L1000 assay39 and to known mutational signatures (see below) that used different indel types. The most drug targets40. favorable indel types according to that benchmark were the simplest where we The drug response data (IC50) for 1502 drugs was downloaded from differentiate 4 indel types: deletions with microhomology (Del-MH), deletions at the PRISM46 database (secondary-screen-dose-response-curve-parameters.csv). microsatellite loci (Del-MS), other deletions (Del-Other), and insertions at any Drugs and cell lines from PRISM and GDSC databases were matched via drug locus (Insertion). names, obtaining in total 348 cell lines and 178 drugs that overlap between the two databases. The gene-level fitness effects for 16,827 genes in 786 cancer cell lines were Simulated cell line exomes. We used TCGA exomes to simulate cell line exomes downloaded as the Integrated cancer dependency dataset from Wellcome Sanger (or more precisely, the variant calls that would originate from a cell line exome) in Institute (release 1) and Broad Institute (19Q3) order to benchmark our ancestry matching procedure for efficacy in reconstructing (“integrated_Sanger_Broad_essentiality_matrices_20200402.zip”) from the the mutation spectrum of cancer cell line exomes. Additionally, we used the Project SCORE database48 and matched to the GDSC cell lines via cell line names, benchmark to optimize the number of clusters in the inference of subpopulations. obtaining in total 517 overlapping cell lines (https://score.depmap.sanger.ac.uk/ To simulate variant calls from cell line exomes, we performed the variant calling for downloads). tumor samples in the same way as for the cell lines, i.e., without matched normal samples. In addition, for every sample we merged SNV calls obtained thus with the somatic SNV calls of tumors, ensuring that the true somatic mutations are a part of Matching of cancer cell lines and primary tumors (ancestry matching pro- a simulated cell line to enable an accurate estimation of errors. In the tumor types cedure). To ensure that genomic data is comparable between cell lines and tumors that were represented both in the set of cancer cell lines and the set of TCGA we performed the following steps: (i) we used only SNVs detected in the regions of tumors, we randomly selected 450 tumor samples taking approximately the same the exome with well sequencing coverage (≥20 reads in 9 out of 10 randomly number of samples per tumor type as in cell lines. In the subsequent analysis selected samples); (ii) to avoid gender bias during analysis, variants on X and Y involving simulated cell line exomes, we removed their corresponding normal chromosomes were not used; (iii) only uniquely mappable regions of the genome samples from the pool of germline samples. We performed the ancestry matching were used, as defined by the Umap k36 criterion;95 (iv) we discarded regions with procedure involving simulated cell line exomes in the same way as described above frequent somatic deletions, namely pan-cancer deletions of significant SCNA96 and (Fig. 1g) and evaluated the accuracy of mutation spectrum reconstruction by frequently deleted chromosome arms (deleted in >18% of tumor samples);97 (v) we comparing the mutation spectrum of a tumor (ground truth) to the corresponding detected copy number changes of cell lines with CNVkit98 and removed deleted simulated cell line mutation spectrum obtained by ancestry matching. As an regions (log2 score < −0.3). From the remaining regions of the exome we selected accuracy measure, we used Absolute Error (the sum of absolute values of differ- common germline variants across the cell line exomes and the TCGA exomes ences between the 96 components of the ground truth spectrum and the recon- (MAF>5% in GnomAD).To perform ‘ancestry matching’ of cancer cell lines with structed spectrum). To estimate the optimal number of subpopulations we varied TCGA normal samples, we performed Principal Component Analysis (PCA) on the number of clusters in the tclust algorithm from 2 to 20, with 13 being the the matrix of common germline variants, followed by clustering according to optimal according to the average Absolute Error of 450 simulated cell lines principal components (PC) (see below). (Fig. 1c). To investigate the variation between germline mutational spectra across human Clustering of cell lines and TCGA germline samples. We employed robust (sub)populations (as reported recently34), we performed hierarchical clustering clustering (“tclust” algorithm;99 discards outlying samples) on the first 140 prin- (based on cosine similarity; stats package in R) of the median mutation spectra of cipal components derived from the common germline variants. Initially, we con- TCGA subpopulations (i.e., median mutation spectra of the TCGA normal samples sidered the first 150 principal components, however, some PCs were attributed to within clusters) (Supplementary Fig. 1E). In addition, we performed PC analysis of the batch effect, i.e., they separate cell lines from TCGA samples. We removed the mutation spectra of TCGA normal samples showing that the main principal top 10 such PCs as determined by feature importance of the random forest clas- component separates the main ethnicity groups (Supplementary Fig. 1D). We note sifier (“randomForestSRC” R package) trained to distinguish between TCGA that, in the previous steps, clustering to infer subpopulations was performed via samples and cell lines. We used the remaining 140 PCs as an input to the tclust common variants (see above), and the trinucleotide spectra were determined after algorithm. Outlying samples, as determined by tclust, were discarded. We varied excluding the non-rare variants (MAF >0.001% in the gnomAD population the number of clusters from 4 to 20, determining the optimal number of clusters database and variants that appear in >5% of samples in our TCGA data set) ruling (13) by using simulated cell line exomes (Fig. 1c; see below). The cancer cell lines out circularity. were matched to their ‘ancestry-matched’ TCGA normal samples, i.e., those Next, we compared our ancestry matching procedure (involving 13 clusters) to belonging to the same cluster as the cell line. We then used ‘ancestors’ to augment the baseline method of determining somatic variants in a cell line genome (i.e., the filtering of germline SNVs from cancer cell line exomes, in addition to filtering filtering of germline variants according to MAF, where SNVs with MAF >0.001% according to the common practices to retain only variants absent or present in very and that appeared in >5% of samples in our dataset are removed) and the bootstrap low frequencies in population databases (see below). self-similarity method (Fig. 1d). The bootstrap self-similarity error was calculated as an error between true and randomly perturbed somatic mutation spectrum, averaged over 100 random runs. For random perturbation, we used the “sampling” Filtering of germline SNVs from cell lines. Here we considered SNVs from the function of the “UPmultinomial” R package. regions of the exome as specified before (steps i–v), with the difference that we used We compared the true number of somatic mutations versus the number SNVs with the sequencing coverage of ≥8 reads in at least 90% samples (9 out of 10 obtained with the ancestry-matching and the MAF filtering (removal of SNVs with randomly selected samples). We next filtered germline variants from cancer cell MAF where MAF>0.001% and that appeared in >5% of samples) to quantify the lines following common practices29,31,79 – here, removing population variants degree of overfiltering and underfiltering of these two filtering approaches (found at MAF>0.001% in the gnomAD population database); additionally, we (Supplementary Fig. 1h). If after filtering, the simulated cell line exome has less filtered germline variants that appeared >5% of samples in our TCGA data set or mutations than the number of true somatic mutations, we consider it overfiltered cell lines data set which removes germline variants that might be particular to this (i.e., some of the true somatic mutations were removed). Otherwise, if it contains data set and also suspected sequencing artifacts). Next, from the remaining SNVs, more mutations than the number of true somatic mutations, we consider it for each cell line and TCGA germline sample, we calculated its trinucleotide underfiltered (i.e., residual germline variants remained). Next, we compared the 14 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE accuracy of mutational signature extraction (see below for the method) with the average per-cancer-type exposure profile (for cancer types available in both cell line trinucleotide profiles of simulated cell lines obtained with MAF filtering and the and PCAWG data). The spectrum-cosine and exposure-cosine similarities were ‘ancestry matching’ approach. We compared the 20 known mutational signatures combined into a single metric given the weight ‘w’ controlling the relative (>0.85 cosine similarity to PCAWG signatures) that were extracted from both sets contribution of each: w×spectrum-cosine + (1–w) × exposure-cosine. We tested of trinucleotide profiles to the mutational signatures extracted from the true different weights: w = 0.5, 0.7, 0.9, and 0.99. For each PCAWG signature, from a somatic trinucleotide profiles. We compared the cosine similarity of trinucleotide set of cell line signatures that match it with spectrum-cosine ≥0.85, we selected the composition of the signatures and the cosine similarity of the exposure scores. The best matching cell line signature as the one with the highest combined metric. ‘ancestry matching’ signatures show statistically significant closer match to the true We considered using additional regression post-processing for assigning somatic signatures than ‘MAF filtering’ signatures in both trinucleotide signature exposures to samples, similarly as done by Petljak et al.31 and Alexandrov composition and exposure scores (Supplementary Fig. 1g). P-values were calculated et al.9. We used the “sigproSS” tool of SigProfilerExtractor framework101 to by Wilcoxon rank-sum test (artifact signatures were excluded from the test). attribute exposures to the extracted cell line signatures (i.e., signatures were used as an input to sigproSS). Additionally, we used a custom script based on regularized regression from “glmnet” R package enforcing different degrees of sparseness. We Extraction of mutational signatures. We extracted cancer cell line mutational considered Ridge, Lasso and Elastic Net regression to model cell line’s mutation signatures from 96 component trinucleotide mutation spectra of cancer cell lines spectra as a combination of cell line signatures, where the resulting regression obtained with the ancestry matching procedure. In total, we used mutation spectra coefficients are considered as exposures to signatures. We enforced non-negative of 930 cancer cell line exomes (Supplementary Data S13). Note that some cell lines coefficients, no interaction term, and used crossvalidation to determine the optimal were excluded because they were assigned to the outlier cluster of the tclust value of the lambda parameter. algorithm, based on their common germline variants (see above). To extract The different procedures for mutational signature extraction we implemented mutational signatures we used a custom R script based on non-negative matrix were evaluated according to the above-described criteria (a–c), results are presented factorization (NMF) as described by Alexandrov et al.33. We additionally imple- in Supplementary Data S1. As a final set of signatures we chose signatures obtained mented a number of signature extraction procedures proposed in the literature with hierarchical extraction (using 0.97 cosine similarity threshold for sample recently9,38,42 and benchmarked the resulting signatures according to the several removal), using 4 indel features in addition to 96 trinucleotide types. We criteria: (a) The agreement with an established set of SBS signatures (PCAWG considered both the cosine similarity in the trinucleotide spectrum and also the signatures9) measured as the number of PCAWG signatures that were recapitu- cosine similarity in the exposures across tissues (in tumor genomes), using w = 0.9 lated at ≥0.85 cosine similarity cutoff; (b) The similarity of the cell line signature as weight on the trinucleotide-cosine and 0.1 on exposure-cosine, for matching our exposure profiles across cancer types to the exposure profiles of PCAWG sig- cell line signatures with the known COSMIC tumor signatures. No post-processing natures, measured as cosine similarity between an average exposure per-cancer- of exposures to signatures was performed since it did not yield improvements type profile of a cell line signature and its matching PCAWG signature (for sig- according to the evaluation criteria, thus, we used the raw NMF scores of the natures recapitulated at ≥0.85 cosine similarity); (c) The accuracy of the signatures exposure matrix. A notable modification from the signature extraction method in predicting drug sensitivity profiles across cancer cell line panels (see below for presented in Alexandrov et al.33 is that we did not limit to a single value of the method). “number of signatures” NMF parameter (choosing based on measures of fit and From the matrix containing mutation spectra of samples, we generated 300 consistency, as in Alexandrov et al.33). Instead, inference was run for many values matrices by bootstrap resampling. One bootstrap sample is obtained by applying of this parameter, and the final set of mutational signatures consists of solutions the sampling function of the “UPmultinomial” R package to each sample’s from different values for the ‘number of signatures’ parameter. spectrum. Next, we applied the NMF algorithm to each of the bootstrap samples This procedure yielded 52 signatures (Supplementary Fig. 2). The cell line (“nmf” function of the “nmfgpu4R” R package to obtain different NMF runs; we signatures are named according to the PCAWG signatures they resemble, e.g., the cell used the Multiplicative update rules algorithm100 with 10,000 as the maximal line signature name SBS4/45L denotes that for PCAWG SBS4 this signature was the number of iterations). For each bootstrap sample, we varied the number of closest match (i.e., SBS4 is the primary signature); 45L denotes that the signature also signatures from 2 to 40. Computations were performed on an NVIDIA GeForce resembled PCAWG signature SBS45 (at cosine similarity ≥0.85). The suffix “L” (for RTX 2080 Ti GPU. “like”) denotes 0.85 ≤ cosine similarity < 0.95 (a somewhat less-close match), while the We implemented the following modifications to the above basic procedure absence of the suffix “L” means cosine similarity ≥0.95. Names of signatures other (each tested independently and some in combinations; Supplementary Data S1) to than the primary signature (if present) are ordered by decreasing cosine similarity. obtain the candidate cell line signatures which were matched to the known tumor We checked if the trinucleotide composition of the exome ‘territory’ covered in PCAWG signatures (see below): (i) Similarly as done in the seminal work of WES sequencing we used could affect the cosine similarities to the known PCAWG Alexandrov et al.33, we clustered each batch of NMF solutions and used the cluster signatures, since they were extracted from the WGS data. We adjusted the medoids as candidate cell line signatures. We used the clara function in clusters R trinucleotide spectra of our signatures to match the WGS spectra and re-calculated package, with Euclidean distance and standardization and “pamLike” options; the cosine similarities to the PCAWG signatures. The cosine similarities are not number of samples to be drawn from the dataset was set to 10% of the number of substantially affected – our signatures map to PCAWG signatures largely samples. Each batch of NMF solutions obtained as described before (one batch irrespective of the adjustment for territories (Supplementary Fig. 17b). Note the corresponds to 300 x n solutions, where n is the number of signatures; varies from one signature that does change between the two normalizations (0.96 vs 0.87) is an 2 to 40) was clustered into k clusters with k-means clustering, varying k from 2 to SBS49-like signature, where SBS49 was previously suggested to be an artifact9. 30. (ii) We applied the RTOL-filtering (relative tolerance) as proposed by Furthermore, to check that the mutational signatures extraction was not biased in Degasperi et al.38, where (separately for each different number of signatures terms of the trinucleotide composition due to the removal of lowly covered regions parameter) NMF runs that diverge for more than 0.1% from the best NMF run are and the common germline loci, we compared the trinucleotide compositions in removed, as measured by the root-mean-square deviation of the factorization. (iii) an exome, our examined territory (largely the exome with lowly-covered regions We implemented the ‘hierarchical extraction’ procedure proposed by Alexandrov removed), and in our examined territory after having removed common SNP loci. et al.9 where NMF is repeated iteratively while removing (or down weighting) the The trinucleotide composition was highly similar (cosine similarities >0.993 in all well-reconstructed samples (cosine similarity above a specified threshold) from a three possible pairwise comparisons; Supplementary Fig. 17a). previous iteration to uncover additional signatures. We allowed a maximum of 3 We next searched for cell line-specific signatures, i.e., signatures that commonly iterations and tested different cosine similarity thresholds ranging from 0.95 to appear in cell line data and do not resemble any of the known tumor signatures. To 0.99. In the case of down weighting, we multiplied the sample’s mutation spectra by this end, we employed k-means clustering (“clara” function in “clusters” R package, 0.05 instead of removing it. (iv) Similarly as in Ghandi et al.42, we performed joint with Euclidean distance, standardization, “pamLike” options and 10% as the number signatures inference of cell line exomes and tumor exomes. To the initial matrix of samples to be drawn from the dataset). Each batch of NMF solutions (from the containing mutation spectra of cell lines, we added mutation spectra of 6154 TCGA signature extraction method selected as final by the evaluation) was clustered into k somatic exomes which were preprocessed in the same way as cell line exomes: the clusters with k-means clustering varying k from 2 to 40. We chose the clustering result same regions of the exome were used (see above for filtering by coverage, (i.e., a set of signatures) where the agreement with PCAWG signatures was mappability, etc.). maximized in terms of the number of cluster medoids that resemble PCAWG From the candidate cancer cell line mutational signatures obtained by the above signatures (at cosine similarity ≥0.85). From such a set of signatures we selected the described methodologies, we searched for the signatures that closely resemble the ones dissimilar from any of PCAWG signatures (cosine similarity <0.8), yielding in ones that were previously found in human cancers9 (referred to as PCAWG total 5 cancer cell line-specific signatures (named SBS-CL). These SBS-CL signatures signatures in the text). We compared the individual signatures to the PCAWG appear together with real signatures and are robust (i.e., they are cluster medoids) signatures searching for the closest matching cell line signature for each PCAWG therefore they are also likely bona fide mutational signatures that might originate, for signature. As a final set of tumor-like cell line signatures, we kept the closest instance, from cell-line specific mutational processes. matching cell line signature if its cosine similarity to the best matching PCAWG In total, this yielded 57 cancer cell line signatures (52 corresponding to known signature exceeded 0.85. For cell line signatures that used indel features in addition tumor signatures, and 5 additional cell line-specific signatures) (Supplementary to 96 tri-nucleotide spectra, cosine similarity was calculated only on 96 spectra Fig. 2). To investigate if some of the extracted signatures are a result of incomplete since the known PCAWG SBS mutational signatures do not have an indel separation of other signatures we considered the coefficients of the Lasso regression component in their current implementation (v3.2). We considered an additional (“glmnet” R package) (Supplementary Fig. 6), where we modeled cell line criterion for matching cell line and tumor signatures where in addition to the signatures extracted in this work as a linear combination of PCAWG signatures mutation-spectrum cosine similarity, a cosine similarity between signature’s (Supplementary Fig. 6). NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 Predicting drug response. We built predictive models of drug response using the carcinoma and lung squamous cell carcinoma (denoted as HNSC/LUSC). We Random Forest (RF) algorithm as implemented in the “randomForestSRC” R additionally considered two groups of cell lines obtained by dividing colorectal cell package. We compared the predictive performance of different predictors: muta- lines according to microsatellite instability status (denoted as COREAD_MSI and tions in cancer driver genes (Muts), recurrent copy number alterations (CNAs), COREAD_MSS; other cancer types did not have enough cell lines with DNA hypermethylation (HypMet), gene expression, previously reported cancer cell microsatellite (in)stability labels to warrant such division). Prior to the association line mutational signatures31,32,42, and mutational signatures extracted in this work. search, we removed 21 cell lines that exhibited either sensitivity or resistance We used Muts, CNAs, and HypMet data as reported by Iorio et al.30, namely, nonspecifically towards a large number of drugs (15 cell lines reported by Abbas- mutations in cancer genes associated with positive selection in tumors, focal Aghababazadeh et al.104 and 6 outlier cell lines considering median ln IC50), as recurrently aberrant copy number segments, and hypermethylated informative 5′ well as additional 29 cell lines that were reported to be misclassified105. In addition, C-phosphate-G-3′ sites in gene promoters. The dependent variable was the con- for each cancer type, we removed outlier cell lines by the total number of mutations tinuous value of a response to a drug (log IC50); we iterated this over all drugs. (remaining after the filtering with the ‘ancestry matching’ procedure; see above), Another possible way to run Random Forest, not employed here, would be to here defined as having the number of mutations >3× interquartile range + upper binarize the log IC50 drug response and run the RF algorithm in classification quartile, or <3 × interquartile range - lower quartile (calculated for each cancer type mode rather than regression mode. However, binarizing the drug response implies separately). We used binarized exposures to mutational signatures where, for each some loss of information, and moreover, it would necessitate an extra parameter in signature, values below the 5% of the value of the second-highest exposure score the analysis (choice of the binarization threshold) so we opted to use continuous (used in order to avoid single high-score outliers observed in some signatures) of drug response values here. that signature were set to 0, and the rest to 1 (Supplementary Data S4). We We built cancer type specific models for each drug separately, considering only empirically tested several different thresholds ranging from 1% to 20% (of the models where at least 15 cell lines with drug response were available to train the second-highest exposure across cell lines), and measured (i) the sparsity of model. For model validation we used 10-fold cross validation which was repeated binarized signatures and (ii) the distribution of binarized exposures across tissues five times to get more stable results. We built random forests with 100 trees and a for signatures with known etiology that should dominantly appear in certain tissues minimal number of samples in terminal nodes set to 2. We assessed the predictive (the “UV” signatures 7a and 7b in skin and the “tobacco” signature 4 in lung). We performance of the predictors by the relative root-mean-square-error (RRMSE), chose 5% since it offered a good tradeoff on these two metrics; sparsity is not too i.e., a root-mean-square error divided by the root-mean-square error of the default high, while the tissue distribution of signatures 7a, 7b, and 4 are reasonable model. The default model predicts a constant value for all cell lines equal to the (Supplementary Fig. 15B). We used normalized (relative) signature exposures, as average ln IC50 across the training set. RRMSE < 1, therefore, denotes better described above. We considered only tests with at least 8 cancer cell lines and at accuracy better than the one of the (uninformative) default model. We define such least 2 non-zero values of a feature. models as predictive. Results are presented as an average RRMSE per cancer type All of the association tests were performed by modeling the drug response (or across all models built for a cancer type (Fig. 2a). Note that, all drugs were not gene fitness score) to associate it with the status of a feature (i.e., a mutational necessarily tested exhaustively across all cancer cell lines, therefore, the number of signature or a genetic feature) searching separately for sensitivity and resistance models per cancer type may differ (one model corresponds to one drug) due to associations. For each association, we calculated the association score as the missing data (<15 cell lines with drug response data available), as well as the minimum (in the case of sensitivity) or maximum (resistance) effect size (Cohen’s number of model per different predictors (mutational signatures, gene expression, d) of the two independent datasets i.e. positive Cohen’s d implies sensitivity and oncogenic mutations, copy number alterations, DNA hypermethylation) due to negative resistance associations. Here, effect size is the Cohen’s d statistic: data availability. In the analysis of complementarity between different predictor difference of mean drug sensitivity (ln IC50) between the cell lines having the types (Fig. 2d, e and Supplementary Fig. 7b) we thus report relative numbers of feature and those not having it, divided by the pooled standard deviation of the predictive models to facilitate fair comparison. data. To obtain the association’s empirical p-value, we performed a randomization To assess the statistical significance of the differences in the predictive procedure where we calculated the association score 100,000 times for the performance (Fig. 2b, c), we follow the recommendations given by Demšar102. randomly shuffled features. For the sensitivity test, the formula for the p value is: More specifically, to statistically compare the predictive performance of multiple p ¼ random score >¼ observed score num:of randomizations and for the resistance is: p ¼ random score <¼ observed score num: of randomizations . predictors over multiple datasets we use the corrected Friedman test and the post- Due to the computational burden, we performed the randomization procedure hoc Nemenyi test103. Here, a dataset corresponds to a pair of one cancer-type and only for associations that had an effect size >0.2 in the primary test. The empirical one drug. Due to the requirements of this test, only the intersection of drugs p-values were adjusted with the Tibshirani-Storey method47. modeled across all cancer types and predictors were considered (due to missing Recent work suggested that MMR-failure signatures can be grouped into a few data the number of drugs modeled per cancer type and/or predictor may differ; see broad types38,106,107, in particular a group enriched with C > T changes, and a group above). For each drug-cancer-type pair, predictors are ranked according to their enriched with T>C (equivalently, A > G) changes. Based on this we considered RRMSE for that drug-cancer-type pair, where rank 1 corresponds to the best (i.e. aggregated mis-match mutational signatures where we merged p-values and effect the lowest) RRMSE, 2 to second-best, etc. Ranks are then averaged across all drug- sizes of SBS6, SBS15, and SBS44L (denoted as SBS-MMR1); SBS21 and SBS26/12L cancer-type pairs to obtain the average rankings of predictors. The Nemenyi test (denoted as SBS-MMR2). Similarly, we considered the aggregate of the two APOBEC performs a pairwise comparison of predictors’ performance based on the absolute signatures SBS2 and SBS13 (denoted as SBS-APOBEC). Pooled p-values of difference of the average rankings of the predictors to assess whether there are aggregated signatures were obtained by Fisher's method, while pooled effect size was statistically significant differences among predictors. The test determines the obtained by averaging. Pooled p-values were adjusted the same as described above. critical difference (CD) for a given significance level α, if the difference between the We consider a two-way association as statistically significant if FDR<15% and average ranks of two predictors is greater than CD, the null hypothesis (that the additionally we imposed an effect size threshold of Cohen’s d > 0.5 in both tests predictors have the same performance) is rejected., i.e., there is a statistically (based on a known set of positive control associations (Supplementary Fig. 16) for significant difference between the two predictors. The results from the Nemenyi GDSC/PRISM and GDSC/PSCORE tests, while for the GDSC/GDSC (same target) post-hoc test are presented with an average ranks diagram102. The average test we required Cohen’s d > 1. In addition, we considered only associations coming predictor rankings are depicted on an axis, in such a manner that the best ranking from cancer types where the inflation factor lambda was below 1.3 (Supplementary algorithms are at the left-most side of the diagram. The algorithms that do not Fig. 8). Note that, in supplementary data also associations with lambda > 1.3 are differ significantly (in performance) for a significance level of 0.05 are connected listed. For some analyses we considered an additional set of associations with an with a bold line, therefore, predictors that are not connected are statistically unadjusted p-value threshold of < 0.005 and the same effect size threshold of significantly different according to the test. Cohen’s d > 0.5. As a rule-of-thumb for interpretation, Cohen’s d = 0.2, 0.5 and 0.8 correspond to small, medium and large effect sizes, respectively108. We used Associations with drug response “QCEWAS” R package to calculate the lambda score to estimate the inflation of Randomization test for associations that replicate in independent datasets. We p-values (calculated separately for sensitivity and resistance associations). performed two-way association testing where we searched for robust associations that replicate in two independent datasets (schematic in Fig. 3a). We considered ‘Golden’ and ‘Silver’ sets of high-priority associations. We collated a list of 3911 three different types of two-way tests where we enforced that, for a given drug, an associations (the ‘silver set’; Supplementary Data S10). To make this list, we con- association between a particular feature (a mutational signature or a genetic fea- sidered all associations tested in the three ‘two-way’ replication tests that pass the ture) and the drug response from the GDSC database is replicated: (1) in the permissive criterion of significance (effect size d > 0.5 in both tests (or d > 1 for the PRISM drug screening database (GDSC/PRISM test) for the same drug, or (2) with GDSC/GDSC two-way test) and a nominal p < 0.005). From these, an association another drug from the GDSC database that shares the same molecular target (between a feature and drug within a certain cancer type) was listed in the silver set (GDSC/GDSC (same target) test), or (3) in the Project SCORE CRISPR/Cas9 if it was confirmed in more than one ‘two-way’ replication tests, or was seen in genetic screen as an association with a protein-coding gene fitness score of one of more than one cancer type. We require that at least one of the supporting asso- the drug’s target proteins (GDSC/PSCORE test). ciations has FDR<25%. Additionally, we collated a ‘golden set’ of high-priority list We considered cancer-type specific associations. We required that in both tests of associations involving mutational signatures and cancer functional events which of a two-way test an association is detected in the same cancer type. We merged we consider to be suitable for follow up work: we require that an association some similar cancer types with a small number of cell lines: esophagus carcinoma involving the same drug is supported in at least three cancer types or is replicated and stomach adenocarcinoma (denoted as ESCA/STAD), glioblastoma and brain in all three ‘two-way’ tests where at least one association has FDR < 25% (995 lower grade glioma (denoted as GBM/LGG), and head and neck squamous cell associations; Supplementary Data S8). 16 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE Abbreviations of cancer types. A list of abbreviations of cancer types used in this 10. Pour, M. V., Lehner, B. & Supek, F. The impact of rare germline variants on study: ALL/CLL, acute/chronic lymphoblastic leukemia; BLCA, bladder urothelial human somatic mutation processes. bioRxix https://doi.org/10.1101/2021.11. carcinoma; BONE, bone cancer other/not classified further; COREAD, colon 14.468508 (2021). adenocarcinoma and rectum adenocarcinoma; COREAD_MSI, microsatellite 11. Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. instable colon and rectum adenocarcinoma; COREAD_MSS, microsatellite stable Nature 500, 415–421 (2013). colon and rectum adenocarcinoma; ESCA, esophageal carcinoma; ESCA_STAD, 12. Shinbrot, E. et al. Exonuclease mutations in DNA polymerase epsilon reveal esophagus carcinoma and stomach adenocarcinoma; EWING, Ewing's sarcoma; replication strand specific mutation patterns and human origins of replication. GBM, glioblastoma multiforme; GBM_LGG, glioblastoma and brain lower grade Genome Res. 24, 1740–1750 (2014). glioma; HNSC, head and neck squamous cell carcinoma; HNSC_LUSC, head and 13. Briggs, S. & Tomlinson, I. Germline and somatic polymerase ϵ and δ neck squamous cell carcinoma and lung squamous cell carcinoma; KIRC, kidney mutations define a new class of hypermutated colorectal and endometrial renal clear cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade cancers. J. Pathol. 230, 148–153 (2013). glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; 14. Grolleman, J. E. et al. Mutational signature analysis reveals NTHL1 deficiency LUSC, lung squamous cell carcinoma; LYMP, lymphoma; MESO, mesothelioma; to cause a multi-tumor phenotype. Cancer Cell 35, 256–266.e5 (2019). MM, multiple myeloma; NB, neuroblastoma; OV, ovarian serous cystadenocarci- 15. Viel, A. et al. A specific mutational signature associated with DNA noma; PAAD, pancreatic adenocarcinoma; SARC, sarcoma other/not classified 8-oxoguanine persistence in MUTYH-defective colorectal cancer. further; SCLC, small cell lung cancer; SKCM, skin cutaneous melanoma; STAD, EBioMedicine 20, 39 (2017). stomach adenocarcinoma; THCA, thyroid carcinoma; UCEC, uterine corpus 16. Campbell, P. J. et al. Pan-cancer analysis of whole genomes. Nature 578, endometrial carcinoma. 82–93 (2020). 17. Kim, J. et al. Somatic ERCC2 mutations are associated with a distinct genomic Reporting summary. Further information on research design is available in the Nature signature in urothelial tumors. Nat. Genet. 48, 600–606 (2016). Research Reporting Summary linked to this article. 18. Zámborszky, J. et al. Loss of BRCA1 or BRCA2 markedly increases the rate of base substitution mutagenesis and has distinct effects on genomic deletions. Oncogene 36, 746–755 (2017). 19. Olivier, M. et al. Modelling mutational landscapes of human cancers in vitro. Data availability Sci. Rep. 4, 4482 (2014). The data sources used for this study are listed below and described further in the 20. Zou, X. et al. Validating the concept of mutational signatures with isogenic cell Methods. The data resulting from our analyses are available in the Supplementary models. Nat. Commun. 9, 1744 (2018). Material, or are otherwise available from the authors upon request. Source data for 21. Szikriszt, B. et al. A comprehensive survey of the mutagenic impact of figures are provided with this paper. Data sources used: WES bam files for human cancer common cancer cytotoxics. Genome Biol. 17, 99 (2016). cell lines (EGA ID number EGAD00001001039, restricted access that can be applied to 22. Macintyre, G. et al. Copy-number signatures and mutational processes in following instructions on EGA), WES bam files for tumors and their matched normals ovarian carcinoma. Nat. Genet. 50, 1262–1270 (2018). (dbGaP accession ID phs000178 [https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/ 23. Polak, P. et al. A mutational signature reveals alterations underlying deficient study.cgi?study_id=phs000178.v11.p8], restricted access that can be applied to following homologous recombination repair in breast cancer. Nat. Genet. 49, 1476–1486 instructions on dbGaP; bams downloaded from NCI Genomic Data Commons [https:// (2017). portal.gdc.cancer.gov/]), drug response data for human cancer cell lines [https://www. 24. Alexandrov, L. B., Nik-Zainal, S., Siu, H. C., Leung, S. Y. & Stratton, M. R. A cancerrxgene.org/]; Release 8.3), PRISM Repurposing dataset 19Q4 [https://depmap.org/ mutational signature in gastric cancer suggests therapeutic strategies. Nat. portal/download/all/], Project Score CRISPR genetic screening data [https://score. Commun. 6, 8683 (2015). depmap.sanger.ac.uk/downloads]; Integrated cancer dependency dataset from Wellcome 25. Jonsson, P. et al. Tumour lineage shapes BRCA-mediated phenotypes. Nature Sanger Institute (release 1) and Broad Institute (19Q3)) [https://score.depmap.sanger.ac. 571, 576–579 (2019). uk/downloads]. Source data are provided with this paper. 26. Gulhan, D. C., Lee, J. J.-K., Melloni, G. E. M., Cortés-Ciriano, I. & Park, P. J. Detecting the mutational signature of homologous recombination deficiency in clinical samples. Nat. Genet. 51, 912–919 (2019). Code availability 27. Davies, H. et al. HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures. Nat. Med. 23, 517–525 (2017). The computer code is available at the GitHub repository: https://github.com/jlevatic/ 28. Garnett, M. J. et al. Systematic identification of genomic markers of drug CellLineMutSigs (https://doi.org/10.5281/zenodo.6463399)109. sensitivity in cancer cells. Nature 483, 570–575 (2012). 29. Barretina, J. et al. The Cancer Cell Line Encyclopedia enables predictive Received: 19 May 2021; Accepted: 9 May 2022; modelling of anticancer drug sensitivity. Nature 483, 603–607 (2012). 30. Iorio, F. et al. A landscape of pharmacogenomic interactions in cancer. Cell 166, 740–754 (2016). 31. Petljak, M. et al. Characterizing mutational signatures in human cancer cancer cell lines reveals episodic APOBEC mutagenesis. Cell 176, 1282–1294.e20 (2019). References 32. Jarvis, M. C., Ebrahimi, D., Temiz, N. A. & Harris, R. S. Mutation signatures 1. Kelley, M. R., Logsdon, D. & Fishel, M. L. Targeting DNA repair pathways for including APOBEC in cancer cell lines. JNCI Cancer Spectr. 2, pky002 (2018). cancer treatment: what’s new? Future Oncol. Lond. Engl. 10, 1215–1237 33. Alexandrov, L. B., Nik-Zainal, S., Wedge, D. C., Campbell, P. J. & Stratton, M. (2014). R. Deciphering signatures of mutational processes operative in human cancer. 2. Gavande, N. S. et al. DNA repair targeted therapy: the past or future of cancer Cell Rep. 3, 246–259 (2013). treatment? Pharmacol. Ther. 160, 65–83 (2016). 34. Harris, K. & Pritchard, J. K. Rapid evolution of the human mutation spectrum. 3. Klinakis, A., Karagiannis, D. & Rampias, T. Targeting DNA repair in cancer: eLife 6, e24284 (2017). current state and novel approaches. Cell. Mol. Life Sci. 77, 677–703 (2020). 35. Mathieson, I. & Reich, D. Differences in the rare variant spectrum among 4. Esteller, M. Epigenetic lesions causing genetic lesions in human cancer: human populations. PLoS Genet. 13, e1006581 (2017). promoter hypermethylation of DNA repair genes. Eur. J. Cancer 36, 36. Aikens, R. C., Johnson, K. E. & Voight, B. F. Signals of variation in human 2294–2300 (2000). mutation rate at multiple levels of sequence context. Mol. Biol. Evol. 36, 5. Lahtz, C. & Pfeifer, G. P. Epigenetic changes of DNA repair genes in cancer. J. 955–965 (2019). Mol. Cell Biol. 3, 51–58 (2011). 37. Yang, W. et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for 6. Knijnenburg, T. A. et al. Genomic and molecular landscape of DNA damage therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, repair deficiency across The Cancer Genome Atlas. Cell Rep. 23, 239–254.e6 D955–D961 (2013). (2018). 38. Degasperi, A. et al. A practical framework and online tool for mutational 7. Diouf, B. et al. Somatic deletions of genes regulating MSH2 protein stability signature analyses show intertissue variation and driver dependencies. Nat. cause DNA mismatch repair deficiency and drug resistance in human Cancer 1, 249–263 (2020). leukemia cells. Nat. Med. 17, 1298–1303 (2011). 39. Subramanian, A. et al. A Next Generation Connectivity Map: L1000 platform 8. Park, S., Supek, F. & Lehner, B. Systematic discovery of germline cancer and the first 1,000,000 profiles. Cell 171, 1437–1452.e17 (2017). predisposition genes through the identification of somatic second hits. Nat. 40. Santos, R. et al. A comprehensive map of molecular drug targets. Nat. Rev. Commun. 9, 2601 (2018). Drug Discov. 16, 19–34 (2017). 9. Alexandrov, L. B. et al. The repertoire of mutational signatures in human 41. Dempster, J. M. et al. Gene expression has more power for predicting in vitro cancer. Nature 578, 94–101 (2020). cancer cell vulnerabilities than genomics. bioRxiv https://doi.org/10.1101/ 2020.02.21.959627 (2020). NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 17 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 42. Ghandi, M. et al. Next-generation characterization of the Cancer Cell Line 74. Lowenthal, B. M. et al. Loss of ARID1A expression is associated with DNA Encyclopedia. Nature 569, 503–508 (2019). mismatch repair protein deficiency and favorable prognosis in advanced stage 43. Haibe-Kains, B. et al. Inconsistency in large pharmacogenomic studies. Nature surgically resected esophageal adenocarcinoma. Hum. Pathol. 94, 1–10 (2019). 504, 389–393 (2013). 75. Ge, H. et al. Mismatch repair deficiency is associated with specific 44. Haverty, P. M. et al. Reproducible pharmacogenomic profiling of cancer cell morphologic features and frequent loss of ARID1A expression in ovarian clear line panels. Nature 533, 333–337 (2016). cell carcinoma. Diagn. Pathol. 16, 12 (2021). 45. Stransky, N. et al. Pharmacogenomic agreement between two cancer cell line 76. Liccardi, G., Hartley, J. A. & Hochhauser, D. EGFR nuclear translocation data sets. Nature 528, 84–87 (2015). modulates DNA repair following cisplatin and ionizing radiation treatment. 46. Yu, C. et al. High-throughput identification of genotype-specific cancer Cancer Res. 71, 1103–1114 (2011). vulnerabilities in mixtures of barcoded tumor cell lines. Nat. Biotechnol. 34, 77. Chou, R.-H. et al. EGFR modulates DNA synthesis and repair through Tyr 419–423 (2016). phosphorylation of histone H4. Dev. Cell 30, 224–237 (2014). 47. Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. 78. Boras, I., Nasser, R., Sabatinos, S. & Antonescu, C. N. Signaling by the Proc. Natl Acad. Sci. USA 100, 9440–9445 (2003). epidermal growth factor receptor regulates DNA repair. FASEB J. 33, 48. Behan, F. M. et al. Prioritization of cancer therapeutic targets using 457.2–457.2 (2019). CRISPR–Cas9 screens. Nature 568, 511–516 (2019). 79. Basu, A. et al. An interactive resource to identify cancer genetic and lineage 49. Lennerz, J. K. et al. MET amplification identifies a small and aggressive dependencies targeted by small molecules. Cell 154, 1151–1161 (2013). subgroup of esophagogastric adenocarcinoma with evidence of responsiveness 80. Tsherniak, A. et al. Defining a cancer dependency map. Cell 170, 564–576.e16 to crizotinib. J. Clin. Oncol. 29, 4803–4810 (2011). (2017). 50. Aparicio, T. et al. The activity of crizotinib in chemo-refractory met-amplified 81. Koh, G., Degasperi, A., Zou, X., Momen, S. & Nik-Zainal, S. Mutational esophageal and gastric adenocarcinomas: results from the AcSé-Crizotinib signatures: emerging concepts, caveats and clinical applications. Nat. Rev. Program. Target. Oncol. 16, 381–388 (2021). Cancer 21, 619–637 (2021). 51. Gonçalves, E. et al. Drug mechanism-of-action discovery through the integration 82. Ben-David, U. et al. Genetic and transcriptional evolution alters cancer cell of pharmacological and CRISPR screens. Mol. Syst. Biol. 16, e9405 (2020). line drug response. Nature 560, 325–330 (2018). 52. Zehir, A. et al. Mutational landscape of metastatic cancer revealed from 83. Liu, Y. et al. Multi-omic measurements of heterogeneity in HeLa cells across prospective clinical sequencing of 10,000 patients. Nat. Med. 23, 703–713 (2017). laboratories. Nat. Biotechnol. 37, 314–322 (2019). 53. Gao, J. et al. Integrative analysis of complex cancer genomics and clinical 84. Park, S. & Lehner, B. Cancer type-dependent genetic interactions between profiles using the cBioPortal. Sci. Signal. 6, pl1 (2013). cancer driver alterations indicate plasticity of epistasis across cell types. Mol. 54. Yang, Y. et al. Loss of ARID1A promotes proliferation, migration and invasion Syst. Biol. 11, 824 (2015). via the Akt signaling pathway in NPC. Cancer Manag. Res. 11, 4931–4946 85. Rauscher, B. et al. Toward an integrated map of genetic interactions in cancer (2019). cells. Mol. Syst. Biol. 14, e7656 (2018). 55. Lee, D., Yu, E. J., Ham, I.-H., Hur, H. & Kim, Y.-S. AKT inhibition is an 86. Christensen, S. et al. 5-Fluorouracil treatment induces characteristic T>G effective treatment strategy in ARID1A-deficient gastric cancer cells. mutations in human cancer. Nat. Commun. 10, 4571 (2019). OncoTargets Ther. 10, 4153–4159 (2017). 87. Shiraishi, Y., Tremmel, G., Miyano, S. & Stephens, M. A simple model-based 56. Samartzis, E. P. et al. Loss of ARID1A expression sensitizes cancer cells to approach to inferring and visualizing cancer mutation signatures. PLoS Genet. PI3K- and AKT-inhibition. Oncotarget 5, 5295–5303 (2014). 11, e1005657 (2015). 57. Bitler, B. G. et al. Synthetic lethality by targeting EZH2 methyltransferase 88. Vöhringer, H., Hoeck, A. V., Cuppen, E. & Gerstung, M. Learning mutational activity in ARID1A-mutated cancers. Nat. Med. 21, 231–238 (2015). signatures and their multidimensional genomic properties with 58. Chan, E. M. et al. WRN helicase is a synthetic lethal target in microsatellite TensorSignatures. Nat. Commun. 12, 3628 (2021). unstable cancers. Nature 568, 551–556 (2019). 89. Cahill, D. P. et al. Loss of the mismatch repair protein MSH6 in human 59. Lieb, S. et al. Werner syndrome helicase is a selective vulnerability of glioblastomas is associated with tumor progression during temozolomide microsatellite instability-high tumor cells. eLife 8, e43333 (2019). treatment. Clin. Cancer Res. 13, 2038–2045 (2007). 60. Kategaya, L., Perumal, S. K., Hager, J. H. & Belmont, L. D. Werner Syndrome 90. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows- Helicase Is Required for the Survival of Cancer Cells with Microsatellite Wheeler transform. Bioinform. Oxf. Engl. 25, 1754–1760 (2009). Instability. iScience 13, 488–497 (2019). 91. Li, H. et al. The sequence alignment/map format and SAMtools. Bioinform. 61. Jager, M. et al. Deficiency of nucleotide excision repair is associated with Oxf. Engl. 25, 2078–2079 (2009). mutational signature observed in cancer. Genome Res. 29, 1067–1077 (2019). 92. Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic 62. Mas-Ponte, D. & Supek, F. DNA mismatch repair promotes APOBEC3- variants. Nat. Methods 15, 591–594 (2018). mediated diffuse hypermutation in human cancers. Nat. Genet. 52, 958–968 93. Picard Tools - By Broad Institute. http://broadinstitute.github.io/picard/. (2020). 94. Karczewski, K. J. et al. The mutational constraint spectrum quantified from 63. Pilati, C. et al. Mutational signature analysis identifies MUTYH deficiency in variation in 141,456 humans. Nature 581, 434–443 (2020). colorectal cancers and adrenocortical carcinomas. J. Pathol. 242, 10–15 95. Karimzadeh, M., Ernst, C., Kundaje, A. & Hoffman, M. M. Umap and Bismap: (2017). quantifying genome and methylome mappability. Nucleic Acids Res. 46, 64. Scarpa, A. et al. Whole-genome landscape of pancreatic neuroendocrine e120–e120 (2018). tumours. Nature 543, 65–71 (2017). 96. Zack, T. I. et al. Pan-cancer patterns of somatic copy number alteration. Nat. 65. Drost, J. et al. Use of CRISPR-modified human stem cell organoids to study Genet. 45, 1134–1140 (2013). the origin of mutational signatures in cancer. Science 358, 234–238 (2017). 97. Taylor, A. M. et al. Genomic and functional approaches to understanding 66. Derrien, A.-C. et al. Germline MBD4 mutations and predisposition to uveal cancer aneuploidy. Cancer Cell 33, 676–689.e3 (2018). melanoma. J. Natl Cancer Inst. 113, 80–87 (2021). 98. Talevich, E., Shain, A. H., Botton, T. & Bastian, B. C. CNVkit: genome-wide 67. Thibodeau, M. L. et al. Base excision repair deficiency signatures implicate copy number detection and visualization from targeted DNA sequencing. germline and somatic MUTYH aberrations in pancreatic ductal PLoS Comput. Biol. 12, e1004873 (2016). adenocarcinoma and breast cancer oncogenesis. Cold Spring Harb. Mol. Case 99. García-Escudero, L. A., Gordaliza, A., Matrán, C. & Mayo-Iscar, A. A general Stud. 5, a003681 (2019). trimming approach to robust cluster Analysis. Ann. Stat. 36, 1324–1345 68. Mimaki, S. et al. Hypermutation and unique mutational signatures of (2008). occupational cholangiocarcinoma in printing workers exposed to haloalkanes. 100. Lee, D. D. & Seung, H. S. Advances in Neural Information Processing Systems Carcinogenesis 37, 817–826 (2016). 13 (eds. Leen, T. K., Dietterich, T. G. & Tresp, V.) 556–562 (MIT Press, 2001). 69. Stiborová, M., Arlt, V. M. & Schmeiser, H. H. DNA adducts formed by 101. Islam, S. M. A. et al. Uncovering novel mutational signatures by de novo aristolochic acid are unique biomarkers of exposure and explain the initiation extraction with SigProfilerExtractor. bioRxiv https://doi.org/10.1101/2020.12. phase of upper urothelial cancer. Int. J. Mol. Sci. 18, 2144 (2017). 13.422570 (2020). 70. Martorana, F. et al. AKT inhibitors: new weapons in the fight against breast 102. Demšar, J. Statistical comparisons of classifiers over multiple data sets. J. cancer? Front. Pharmacol. 12, 662232 (2021). Mach. Learn. Res. 7, 1–30 (2006). 71. Song, M., Bode, A. M., Dong, Z. & Lee, M.-H. AKT as a therapeutic target for 103. Nemenyi, P. Distribution-free Multiple Comparisons (Princeton University, cancer. Cancer Res. 79, 1019–1031 (2019). 1963). 72. Shen, J. et al. ARID1A deficiency promotes mutability and potentiates 104. Abbas-Aghababazadeh, F., Lu, P. & Fridley, B. L. Nonlinear mixed-effects therapeutic antitumor immunity unleashed by immune checkpoint blockade. models for modeling in vitro drug response data to determine problematic Nat. Med. 24, 556–562 (2018). cancer cell lines. Sci. Rep. 9, 14421 (2019). 73. Allo, G. et al. ARID1A loss correlates with mismatch repair deficiency and 105. Salvadores, M., Fuster-Tormo, F. & Supek, F. Matching cell lines with cancer intact p53 expression in high-grade endometrial carcinomas. Mod. Pathol. 27, type and subtype of origin via mutational, epigenomic, and transcriptomic 255–261 (2014). patterns. Sci. Adv. 6, eaba1862 (2020). 18 NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30582-3 ARTICLE 106. Fang, H., Zhu, X., Oh, J., Barbour, J. A. & Wong, J. W. H. Deficiency in DNA Competing interests mismatch repair of methylation damage is a major mutational process in The authors declare no competing interests. cancer. bioRxiv https://doi.org/10.1101/2020.11.18.388108 (2020). 107. Németh, E. et al. Two main mutational processes operate in the absence of DNA mismatch repair. DNA Repair 89, 102827 (2020). Additional information Supplementary information The online version contains supplementary material 108. Cohen, J. Statistical Power Analysis for the Behavioral Sciences. (Academic available at https://doi.org/10.1038/s41467-022-30582-3. Press, 2013). 109. Levatic, J., Salvadores, M., Fuster-Tormo, F. & Supek, F. Code Related To The Correspondence and requests for materials should be addressed to Fran Supek. Article Mutational Signatures Are Markers Of Drug Sensitivity Of Cancer Cells (Zenodo, 2022). Peer review information Nature Communications thanks Federico Abascal and the 110. Rahbari, R. et al. Timing, rates and spectra of human germline mutation. Nat. other anonymous reviewer(s) for their contribution to the peer review of this work. Genet. 48, 126–133 (2016). Reprints and permission information is available at http://www.nature.com/reprints Acknowledgements Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Work in the F.S. lab was supported by the ERC Starting Grant “HYPER-INSIGHT” published maps and institutional affiliations. (grant number 757700 to F.S.); the Horizon2020 RIA project “DECIDER” (grant number 965193 to F.S.); the Spanish Ministry of Science, Education and Universities project “REPAIRSCAPE” (grant number PID2020-118795GB-I00 to F.S.); the State Agency for Research of the Ministry of Science and Innovation – Severo Ochoa Centre of Excellence Open Access This article is licensed under a Creative Commons Award (grant number CEX2019-000913-S to IRB Barcelona); and the CERCA Gen- Attribution 4.0 International License, which permits use, sharing, eralitat de Catalunya funds (to IRB Barcelona). J.L. was co-funded from the European adaptation, distribution and reproduction in any medium or format, as long as you give Union’s Framework Programme for Research and Innovation “Horizon 2020”, under the appropriate credit to the original author(s) and the source, provide a link to the Creative Marie Skłodowska-Curie PROBIST grant agreement No. 754510 (to J.L.). M.S. was Commons license, and indicate if changes were made. The images or other third party funded by the FPU fellowship of the Spanish government (FPU2017, to M.S.). material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Author contributions regulation or exceeds the permitted use, you will need to obtain permission directly from JL collected and curated the data, devised methods and implemented software, performed the copyright holder. To view a copy of this license, visit http://creativecommons.org/ most analyses, visualized the data, and interpreted the results. MS devised and imple- licenses/by/4.0/. mented the replicated-association tests, as well as performed corresponding ana- lyses. FFT processed cell line genomic data. FS conceptualized the analysis, devised the overall methodology, interpreted results, and supervised the study. FS, JL, and MS drafted © The Author(s) 2022 the manuscript jointly. NATURE COMMUNICATIONS | (2022)13:2926 | https://doi.org/10.1038/s41467-022-30582-3 | www.nature.com/naturecommunications 19