A practical guide for mutational signature analysis in hematological malignancies

Analysis of mutational signatures is becoming routine in cancer genomics, with implications for pathogenesis, classification, prognosis, and even treatment decisions. However, the field lacks a consensus on analysis and result interpretation. Using whole-genome sequencing of multiple myeloma (MM), chronic lymphocytic leukemia (CLL) and acute myeloid leukemia, we compare the performance of public signature analysis tools. We describe caveats and pitfalls of de novo signature extraction and fitting approaches, reporting on common inaccuracies: erroneous signature assignment, identification of localized hyper-mutational processes, overcalling of signatures. We provide reproducible solutions to solve these issues and use orthogonal approaches to validate our results. We show how a comprehensive mutational signature analysis may provide relevant biological insights, reporting evidence of c-AID activity among unmutated CLL cases or the absence of BRCA1/BRCA2-mediated homologous recombination deficiency in a MM cohort. Finally, we propose a general analysis framework to ensure production of accurate and reproducible mutational signature data.

T he advent of next generation sequencing has profoundly changed both the research and clinical approach to cancer in the last 10 years 1 . While the cancer genome landscape may be composed of thousands of events, only a minimal fraction of them can be considered as drivers [2][3][4][5] . Despite the majority of tumor mutations do not have a functional role, the entire coding and non-coding mutational catalog can be extremely informative for the identification of the mutational processes operative in different cancer types during initiation and progression 4,[6][7][8][9][10] .
Historically, a simple analysis of single-nucleotide variants (SNVs) as a six-class mutational spectrum (C•G → A•T, C•G → G•C, C•G → T•A, T•A → A•T, T•A → C•G, and T•A → G•C) has highlighted how different cancer types are characterized by different contributions from each class, some of which strongly associated with distinct exogenous carcinogens exposure 11,12 . For example, the C•G → A•T transversion is related to smoking in lung cancer samples 13 , and the C•G → T•A transition is significantly over-represented in skin cancers related to UV light exposure 11,12,14 . Following on from these preliminary observations, different approaches have been suggested to gain resolution in the analysis of these so called mutational signatures. Combining the six possible SNV classes together with their trinucleotide contexts (i.e., the bases 5ʹ and 3ʹ of the mutated nucleotide) all SNVs have been classified into 96 possible combinations 6,7,15 . This classification has then been used to extract >30 different mutational signatures with a non-negative matrix factorization (NNMF) approach from a large series of whole-genome (WGS) and exome (WES) sequencing data 6,16,17 . Some of these signatures are specifically associated with defects of DNA repair mechanisms, exposure to exogenous carcinogens, or different patterns of structural variants (SVs), suggesting they truly reflect known and unknown mutational processes shaping the genome of each cancer type 10,15,[17][18][19][20] . Further to corroborating their biological relevance, some mutational signatures are also associated with a distinct clinical outcome and emerged as potential biomarkers for novel target therapies 18,19,21,22 .
Since this initial effort, several alternative approaches to NNMF have been proposed to improve the mathematical efficacy and biological accuracy of mutational signatures extraction from the 96-class profile of each cancer 6,7,10,[23][24][25][26][27][28][29] . However, the field of mutational signature extraction still lacks a unanimous consensus and standardization of analysis, often resulting in discrepancies between results from similar datasets obtained using different methodological approaches 4,9,10,21,22,[30][31][32][33] . As WGS and WES are becoming common practice, with implications for both basic and translational research, we believe that more should be done to improve the performance and the reproducibility of mutational signature analysis.
In this study, we use different publicly available bioinformatics tools to analyze public datasets from multiple myeloma (MM) and chronic lymphocytic leukemia (CLL) samples, and validate our findings in additional published and unpublished sequencing data from acute myeloid leukemia (AML) samples, to summarize the main factors that should be considered in a high-confidence mutational signature analysis. We discuss sources of bias and pitfalls, and provide a rational and practical approach that could be validated in other independent studies.

Results
Common issues of mutational signature analysis. All different mutational signature analysis algorithms produce a decomposition matrix C ≈ SE, where C is the catalog matrix, with mutation types as rows and samples as columns, S is the signature matrix, with mutation types as rows and signatures as columns, and E is the exposure matrix, with signatures as rows and samples as columns ( Supplementary Fig. 1). Nevertheless, different approaches can be divided in two main groups: (i) the ones that allow de novo signature extraction (e.g., the NNMF framework from Alexandrov et al.) 6 , where given a matrix C the algorithm finds matrices S and E such that C ≈ SE, and (ii) the ones that fit the 96mutational catalog to a pre-selected list of signatures (e.g., the 30 COSMIC signatures), where given C and S the algorithm finds E such that C ≈ SE. An example of algorithm of the second group is deconstructSigs 24 . Both approaches can be extremely informative in different settings, though it is not always easy to determine when and how to use one or the other. Working with mutational signatures analysis with either group of algorithms, we identified three main issues. The first is the ambiguous signature assignment that occurs when different combinations of signatures can explain equally well the same mutational catalog. This issue may arise when multiple so called flat mutational signatures are potentially present in the same data set (e.g., COSMIC signatures 3, 5, and 8) ( Supplementary Fig. 2) 6,31,34 . The second usually occurrs when localized mutational processes are not investigated. In fact, when a signature extraction is performed using all the mutations found in a genome (or exome), only mutational signatures induced by mutational processes that act across the entire genome are usually identified. Localized mutational processes are often responsible for a small proportion of the total number of genome-wide mutations, and thus are generally missed 9,10,35,36 . The third common issue is the bleeding of signatures. It is biologically sound to assume that each cancer sample presents the activity of a limited number of mutational processes. If an extraction is performed on a heterogeneous set of samples, it is possible that signatures present in only part of the set are also erroneously assigned to the entire set. This is mostly due to the algorithms' assumption that all analyzed samples share a similar mutational signature landscape and to the fact that some signatures are similar to each other.
Mutational signature extraction vs. fitting. As mentioned above, a signature analysis can be performed using either a de novo extraction or a fitting approach based on a pre-selected reference list of known signatures (e.g., the 30 COSMIC signatures).
The first approach extracts recurrent patterns of variants in their trinucleotide context from the input data allowing the unbiased identification of both known and novel mutational processes. However, the weakness of this approach is that extracted signatures often do not appear identical to the reference ones. Common problems are: (i) union of co-occurrent multiple signatures into one; (ii) over splitting of one mutational signature into two or more. All these factors can significantly impact the assignment of extracted signatures to the reference ones 6,31 , and this may introduce bias in the estimation of each signature's activity in the samples.
The second approach fits the input data to a suitable reference list of mutational signatures, allowing a better estimation of each signature's relative and absolute contribution for each sample. However, a fitting approach is not able to discover any novel signature and thus needs a priori knowledge of which mutational processes may be operative in that sample cohort. Furthermore, these approaches may be prone to overfitting leading to signature bleeding, i.e., they may assign all signatures from the reference list to all samples. Therefore, before running any fitting algorithm, it is crucial to have at least some knowledge about which mutational processes are operative in the samples to avoid both false positives (overfitting of signatures) and false negatives (missing novel mutational process).
To provide an example of the problems that a fitting algorithm can pose to the interpretation of data if analyzed without any a priori knowledge, we used a cohort of 30 MM cases (Supplementary Table 1), which have been extensively characterized from a genomic point of view. Here, we first applied NNMF-based, de novo extraction algorithms, i.e., the framework from Alexandrov et al. 6,7 . (Fig. 1a, b) and the NNMF approach of the mutationalPatterns R package 37 (Supplementary Sofware 1). Both NNMF approaches extracted five signatures: the clock-like signatures (Signature 1 and 5 merged together), APOBEC (Signature 2), Signature 8, Signature 9, and a new signature named MM1, again highlighting the impact that NNMF approaches can have in new signature discovery (Supplementary Data 1) 6,9,16,23 . Then, using the same input data we then ran two fitting approaches (deconstructSigs and the fitting approach of mutationalPatterns) without a priori knowledge of the active mutational processes in MM and therefore including all 30 COSMIC signatures. DeconstructSigs forced the extraction of a large number of signatures, including ones not previously extracted by NNMF, and some of which clearly representing false positives ( Fig. 1c and Supplementary Sofware 1). For example, the contribution of tobacco-smoking (COSMIC Signature 4) to MM development can most likely be ruled out, as can the contribution of the liver-specific Signature 16 ( Fig. 1c) 17,31,38 . Furthermore, the new signature MM1 was not identified, simply because it was not included in the COSMIC catalog. To reduce false positives, some corrections can be applied to the fitting approach. For example, deconstructSigs uses forward selection to estimate a minimal number of signatures, and removes a signature's contribution to a sample if it accounts for <6% of the sample's mutations. In contrast, mutationalPatterns fitting approach does not introduce any correction while attempting to fit all 30 COSMIC signatures. In this case, a false-positive minimal contribution of unlikely signatures was detected in all patients (mismatch repair, UV light, tobacco-smoking etc.) (Supplementary Sofware 1). Altogether, this shows that fitting approaches may crucially alter the inferred mutational signature landscape in MM. Conversely, when we ran deconstructSigs and mutationalPatterns imputing the shortlist of COSMIC signatures previously identified by the extraction approaches (i.e., NNMF), this led to a more biologically sound assignment and quantification of the absolute and relative contribution of each process (including the new signature MM1) for each sample, significantly reducing the false-positive signatures ( Fig. 1d and Supplementary Sofware 1).

Absence of BRCA-mediated Homologous Recombination
Deficiency in MM. The genomic profile of MM is characterized by several recurrent and private cytogenetic aberrations, making it one of the most complex hematological malignancies from this point of view 3,21,[39][40][41][42][43][44] . Recently, using a fitting approach like deconstructSigs with default parameters 24 , a potential activity from Signature 3 has been proposed in a significant fraction of MMs 32 . This mutational signature is well-known to correlate with BRCA1 and BRCA2 bi-allelic loss and homologous repair deficiency (HRD) in different solid cancers 6,18,20,45 . Signature 3 was indeed observed in our MMs when either mutationalPatterns or deconstructSig fitting approaches were run using all 30 COMISC signatures ( Fig. 1c and Supplementary Sofware 1), but not observed in our signature extraction.
To positively confirm whether or not signature 3 is present in our samples, we used two validation strategies: (1) determine whether the pattern of Signature 3 is necessary to explain the  mutational patterns observed in the samples; (2) analyze additional genomic features to determine the presence of HRD. First, to establish whether Signature 3 is required to explain the catalog of mutational signatures in our samples, we determined whether including or not Signature 3 in our analysis would affect the reconstruction error, i.e., the difference between the original catalogs and the fitted linear combination of signatures for each sample (see Methods). The inclusion of Signature 3 produced a statistically significant lower reconstruction error (measured as KL divergence, root mean squared error (RMSE) or cosine similarities), which can be attributed to the inclusion of an additional signature in the linear combination. However, the reconstruction error is not qualitatively different in the absence of Signature 3 ( Supplementary Fig. 3a-c, g-i). In contrast, when Signature 3 is used in place of either Signature 8 or Signature 5, we have a qualitative increase in the reconstruction error ( Supplementary Fig. 3d-f, j-l). Interestingly, when Signature 3 is excluded, the mutations that were assigned to Signature 3 seem to be reassigned mostly to the other flat Signatures 8 and 5 ( Supplementary Fig. 4). This evidence indicates that Signature 3 is not necessary to explain the patterns of SNV mutations in the samples. Conversely, Signature 8 and Signature 5 emerged as the most significant processes, and the ones that are likely active.
Next, we used an orthogonal approach to detect the presence of BRCA1/BRCA2-like HRD in our MM samples ( Fig. 2): to this end, we applied the recently published HRDetect tool 18 , a highly accurate classifier that estimates the presence of BRCA1/BRCA2like HRD in solid cancers, trained on multiple mutational patterns, including COSMIC Signature 3, COSMIC Signature 8, microhomology-mediated deletions, Rearrangement Signatures 3 and 5 (unclustered short tandem duplications and deletions, respectively) 20 and the HRD index 46 . If we exclude Signature 3 from our analysis, none of the 30 MM samples would be classified as HRD, as they do not appear to be enriched with the patterns that are typical of the BRCA1/BRCA2-type of HRD: there is a low proportion of microhomology-mediated type of small deletions, the HRD-LOH index 46 (Fig. 2b). This sample, characterized by multiple complex events and chromothripsis 47 , is likely to be a false positive generated by the erroneous inclusion of Signature 3 in our analysis. In fact, it lacked the characteristic unclustered genome-wide rearrangements and predominance of microhomology-mediated type of small deletions (Fig. 3a, b and Supplementary Figs. 5 and 6). Finally, if we included Signature 3, we would expect some correlation between the HRDetect score and the assignment of Signature 3, since they both correlate with HRD. However, such correlation is absent in our analysis (Fig. 2b, c).
In conclusion, fitting approaches like deconstructSigs (or mutational pattern) tend to force the assignment of flat signatures, such as Signature 3, to samples when all 30 COSMIC signatures are used as input (Fig. 1c, Fig. 3a, and Supplementary Sofware 1). However, we demonstrated that Signature 3 is not necessary to explain the mutational patterns of MM samples, which furthermore do not show a genomic landscape consistent with BRCA1/BRCA2 loss and its related HRD in terms of 96-class profiles, number of microhomology-mediated deletions and internal tandem duplications as compared to breast cancer (Fig. 3b, c and Supplementary Figs. 5 and 6). We therefore suggest that Signature 3 (and consequently BRCA1/2-mediated HRD) is not biologically active in our MM samples, and it likely represents a false-positive call. Rather, we believe that the right signatures to be annotated in these samples are Signature 8, widely involved in solid and hematological cancers with an unknown etiology 6 , and Signature 5, a flat clock-like process present in normal and cancer tissues 16 . This of course does not exclude the possibility that a larger cohort of MM samples may show cases of BRCA1/2-like HRD, though again, we have no evidence that this is the case in our cohort.
Localized hypermutation. When a naive B-cell passes through the germinal center (GC), it is usually exposed to the activity of activation-induced cytidine deaminase (AID), which is responsible for a very unique genetic process called somatic hypermutation (SHM) of the B-cell receptor (BCR) variable region (VDJ) 48 . This mutational process plays a critical role in the antibody diversification promoting mutations and aminoacidic changes on immunoglobulin heavy and light chain (IGH/IGK/ IGL) genes in order to increase the B-cell receptor (BCR) affinity to distinct antigens 48 . Chronic lymphocytic leukemia (CLL) is well-known to be characterized by two main biological subgroups: one dependent on GC exposure and one independent (Supplementary Data 2). These are differentially diagnosed by recognizing patterns of AID-driven somatic hypermutation in one group (mutated CLL, M-CLL) and not in the other (unmutated CLL, U-CLL) 5,49-53 . MM and M-CLL are post-GC lymphoproliferative malignancies, and their (pre)malignant cells are exposed to AID activity 9,32 . This mutational process, named canonical-AID (c-AID), has been known for years and is specifically active on IGH/IGK/IGL loci 48,54,55 ; however, thanks to mutational signatures analysis, an alternative AID-driven mutational process has been recently observed genome-wide in all post-GC lymphoproliferative disorders 6,10,52,53 . This process was named non-canonical AID (nc-AID; COSMIC Signature 9) and differs from the above-mentioned c-AID in terms of preferential trinucleotide context, genomic distribution and associated cell cycle phase ( Supplementary Fig. 7) 55 . In contrast to nc-AID, the c-AID signature is generally not identified by de novo signature extraction algorithms because it is localized and its limited activity is diluted below the threshold of detection by the larger number of genome-wide mutations generated by other processes (see the lack of its detection in all MM and CLL samples in Fig. 1, Supplementary Data 1, and Supplementary Sofware 1 and 2) 9,10,52 . However, identification of the mutational burden of c-AID and its aberrant targets (e.g., BCL6 54 ) can be extremely informative to compare the genomic landscape of different lymphoproliferative disorders and their different biological origins. The characterization of this localized mutational process can be performed in two ways, with either extraction or fitting algorithms after inclusion of the c-AID 96-class profile ( Supplementary  Fig. 7), currently not part of the COSMIC panel: (1) Considering only hypermutated regions, i.e., those with >5 mutations with a median inter-mutational distance of < 1 Kb; 6,9,15,47 (2) Considering only mutations that occur within known c-AID targets, in particular the IGH/IGK/IGL loci 52 . Both approaches can identify c-AID in both MMs and CLLs (Fig. 4), i.e., two neoplasms where activity of this enzyme is expected. Interestingly, and confirming other previous preliminary data 10 , c-AID activity was also detected in a fraction of U-CLL patients despite the GCindependent pathogenesis. Specifically, in MM and to a greater extent in M-CLL, >10% of these mutations were observed within coding genes, in particular across the VDJ region of the IGH locus; conversely, among U-CLL this activity involved mostly the non-coding part of the IGH locus, in particular within the class switch recombination loci ( Supplementary Fig. 8a-d). These data are in line with the ability of WES to identify c-AID signature within the IG loci only among M-CLL cases 52 , and strengthen the need for WGS for a comprehensive signature analysis.  18 . Scores are ordered from highest to lowest and a classification threshold of 0.7 is used to classify samples as HRD-positive (see Davies et al. 18 ). Below each score, the contribution of the six features that are used by HRDetect is shown. Each contribution is given by the amount of a feature in a sample, log-transformed and standardized according to mean and standard deviation of the features in Davies et al. 18  Furthermore, in contrast to MM and M-CLL cases, nc-AID was not active in IGH regions from U-CLL cases (Fig. 4). Confirming previous reports on a potential ongoing AID activity in U-CLLs 10 , a significant higher fraction of subclonal c-AID mutations (i.e., late mutations) was observed among this group of CLLs (Supplementary Fig. 8e). Conversely, c-AID mutations were mostly detected at clonal level (i.e., early mutations) in M-CLL and MM, confirming the recently reported decreased AID activity in late stages of these diseases 9,10 . Overall, these data suggest a possible non-VDJ and GC-independent role of c-AID among U-CLLs (Fig. 4) 10,56 .
To better characterize the c-AID activity on known loci, we usually prefer to focus on mutations within known c-AID targets rather than to identify hypermutated regions. In fact, most of c-AID mutations occurred close to different VDJ breakpoints, where distant genomic regions are joined by the RAG/AID complex during early stage of B-cell development before the GC exposure 48 .
This means that inter-mutational genomic distance does not reflect the true position of these mutations and should be corrected for the VDJ structure to identify mutations caused by c-AID activity ( Supplementary Fig. 9). This also applies to localized hypermutation events (i.e., kataegis) around complex structural variants (i.e., chromothripsis), where the cancer chromosomal structure significantly differs from the reference 15,47 .
As mentioned above, this kind of analysis can be also directed on known c-AID aberrant targets, such as BCL6, allowing the characterization of clustered mutational processes active around these critical oncogenes and key GC regulators ( Supplementary  Fig. 10) 54 . In our series, BCL6 was involved in localized mutational processes in M-CLL and MM reflecting their GC exposure, as expected; conversely, U-CLLs did not show any evidence of this process, confirming the GC-independent pathogenesis and suggesting the existence of a GC-unrelated AID activity in this group of patients.    T>A  T>C  T>G   C>A  C>G  C>T  T>A  T>C  T>G   C>A  C>G  C>T  T>A  T>C  T>G   C>A  C>G  C>T  T>A  T>C  T>G   C>A  C>G  C>T  T>A  T>C  T>G   C>A  C>G  C>T  T>A  T>C    where the color of each dot represents the mutation class), indels (dark green = insertion; and brown = deletion); copy number variants (red = deletions, green = gain), rearrangements (blue = inversion, red = deletions, green = ITD, black = translocations). PD26419a is the only patient with a slightly high HRDetect score when analyzed including Signature 3. c Circos plots of a breast cancer sample without BRCA deficiency (PD4069a), one with BRCA1 deficiency (PD6413a) and one with BRCA2 deficiency (PD4954a). The MM genomic landscape shows significant differences to the two BRCA-deficient breast cancers, in particular in terms of numbers of indels and SVs, suggesting BRCA-driven HRD is not present in the MM samples analyzed SHM is only present in post-GC B-cells, however it is not the only example of localized hypermutation in cancer. An instance of localized hypermutation termed kataegis has been found across many cancer types and is often promoted by aberrant activity of the APOBEC family of DNA deaminases 47,57 . We have previously reported widespread and localized activity of APOBEC in MM (Fig. 5a-c) 9 where it is recurrently associated with complex rearrangements such as chromothripsis, similarly to what has been reported in several other solid cancers 47 . Furthermore, here we report the first case of APOBEC-mediated kataegis in a therapy-related AML case, again associated with a complex rearrangement (Fig. 5d-f). Previously, APOBEC was never reported as active in AML 6,31 . Overall, our findings stress the importance of performing ad-hoc signature analysis in localized mutational events, since this can highlight specific pathogenetic mechanisms across different cancer types.
Inter-sample bleeding. Both WGS and WES data have clearly shown that M-CLL samples are characterized by a very distinct mutational process (COSMIC Signature 9), reflective of the genome-wide nc-AID activity within the GC 6,10,52 . Conversely, we would expect the absence of nc-AID signature in U-CLL, as these cases do not develop through the GC. To validate this assumption, we performed a de novo signature extraction on all CLLs, using either the Alexandrov et al. 6 framework or the mutationalPatterns 37 NNMF function (Supplementary Data 1). A nc-AID signature was assigned to all samples, with high activity in M-CLL samples and a much lower contribution in U-CLLs ( Fig. 6 and Supplementary Sofware 2). This represents a typical example of inter-sample bleeding effect caused by the assumption that all these samples shared a similar mutational landscape. This incorrect assignment would not be readily highlighted if the biology underlying CLL pathogenesis was not thoroughly known. To obviate this problem, we propose two approaches. In the first, we re-fit the extracted signatures. Here, signatures are first extracted with a de novo approach. Then, a fitting algorithm such as deconstructSigs is applied using only the signatures extracted by NNMF to clean up low-contribution signatures, mostly representing false positives (Fig. 6b, c). The second approach involves performing separate extractions. NNMF is run independently on two sets of samples, split using prior knowledge of the IGHV mutational status evaluated, for example, by Sanger sequencing (Fig. 6d, e and Supplementary Data 2). Either approach successfully removed the nc-AID signature from U-CLL samples, in accordance with the pathogenesis of this CLL subgroup known not to be exposed to GC activity (Fig. 6d, e) 58 .
This kind of a priori biological and clinical knowledge is not available for all cancer types. However, a simple clustering analysis based on the relative contribution of NNMF-extracted mutational signatures may also highlight the heterogeneity in signature activity and therefore help in the identification of distinct groups of patients, based on exposure to different mutational processes (Supplementary Fig. 11). Next, either a second NNMF run or a fitting approach using the NNMF shortlist can be performed on each single subgroup, as explained above 21 .
This inter-sample bleeding of signatures is of course a universal phenomenon and as such can be also observed in non-B-cell hematological malignancies. To extend the validity of our findings we therefore focused on acute myeloid leukemia (AML), where we (i) performed WGS on two cases of therapyrelated AMLs (t-AML) arisen after platinum-based chemotherapy for ovarian carcinoma and (ii) analyzed publicly available WGS data from the TCGA repository of primary AML cases (n = 50) 59 . In this setting, we extracted four main mutational processes: Signature 1, Signature 5 and two signatures currently not included in COSMIC. Of these, one was recently associated with platinum exposure (platinum signature) and the second to the The platinum signature contributed for >30% of the mutational burden of t-AMLs, but its activity was also found among primary AML from TCGA (Fig. 7c). This is inconsistent with the prior knowledge of these samples being treatment-naive. Confirming that platinum signature in primary AML samples represents a further example of inter-sample bleeding, analysis of TCGA primary AMLs without the two t-AML cases led to disappearance of the Platinum Signature ( Fig. 7d and Supplementary Sofware 3). Furthermore, our analysis confirmed the added benefit of performing a de novo signature extraction as a first approach, as two out of four mutational signatures extracted in this cohort of 52 AMLs are not currently included in COSMIC.

Discussion
In this study, we explored caveats and pitfalls of mutational signature analysis using whole-genome sequencing data from three common hematological neoplasms, focusing on the sample set preparation and post-algorithm interpretation processes. Furthermore, we showed how a comprehensive and detailed mutational signature analysis can provide relevant biological insights within different and well characterized cancer types, such as the c-AID activity among UM-IGHV, the absence of BRCA1/ BRCA2-mediated HRD in a MM cohort and two mutational processes in AML, one related to platinum and one less characterized related to stem and progenitor bone marrow cells 31,38,[60][61][62] .
With the rapid increase in the number of tumor genomes sequenced, novel mutational signatures can be identified using several approaches discussed in this work. However, blind trust on out-of-the-box results from public tools can produce an incomplete representation of signatures, or the inclusion of false positives. Our results contain useful practical considerations that can resolve some of the uncertainty in the use of different algorithms, and in the interpretation of the results. Important caveats and pitfalls a scientist can face in mutational signature analysis can usually be recognized and corrected by a priori knowledge of the biology of the tumor and by deep understanding of the way each algorithm works. For example, in CLL it is known that nc-AID exposure within the germinal center is only present among M-CLL cases. Therefore, the finding of Signature 9 activity in U-CLL must be regarded to as artefactual, related to the bleeding phenomenon that is common among de novo NNMF-based approaches. Knowing weaknesses and strengths of each approach, we proposed solutions to improve the accuracy of signature identification, with results that are biologically plausible. The main point of this study is in fact to highlight how the statistical and mathematical methods are important, but they must be used with expertize and combined with a good knowledge of the cancer type being studied. This is especially true when it comes to assignment of flat signatures: our original analysis demonstrates that the previously identified presence of BRCA1/BRCA2-like HRD in an MM cohort is likely to be a false-positive call of fitting algorithms 32 , but this can only be demonstrated knowing the actual genomic consequences of BRCA deficiency in cancers and comparing them to what is seen in MM. Of course, our results only argue against the presence of BRCA1/BRCA2-type of HRD in our MM cohort, as we and others have convincingly demonstrated that a subset of MM patients are characterized by a significant grade of genomic instability 3,21,22,44,63-65 .      CLL5  CLL12  CLL16  CLL20  CLL23  CLL25  CLL26  CLL27  CLL30  CLL48  CLL56  CLL58  CLL63  CLL82  CLL83  CLL134  CLL138  CLL139  CLL141  CLL145  CLL157  CLL166  CLL176  CLL184  CLL188  CLL277  CLL282  CLL290  CLL294  CLL296  CLL306  CLL318  CLL356  CLL371  CLL386  CLL435  CLL442  CLL723  CLL749  CLL761  CLL853  CLL875  CLL1059  CLL1076  CLL1101  CLL1169  CLL1183  CLL1192  CLL1297  CLL1344  CLL1371  CLL1384  CLL1431  CLL148  CLL802  CLL1432  CLL10  CLL1191  CLL1430  CLL177  CLL358  CLL577  CLL1267  CLL15  CLL6  CLL677  CLL278  CLL128  CLL44  CLL753  CLL1078  CLL9  CLL199  CLL618  CLL880  CLL1446  CLL174  CLL372  CLL661  CLL1319  CLL29  CLL283  CLL1164  CLL367  CLL795  CLL804  CLL1424  CLL192  CLL467  CLL1179  CLL564  CLL745  CLL343  CLL393  CLL39  CLL832  CLL519  CLL1525  CLL783  CLL1360  CLL137  CLL90  CLL342  CLL33  CLL4  CLL654  CLL1353  CLL288  CLL1426  CLL1425  CLL1427  CLL84  CLL1276  CLL151  CLL308  CLL776  CLL244  CLL477  CLL32  CLL862  CLL1056  CLL628  CLL1423  CLL179  CLL473  CLL594  CLL523  CLL684  CLL1477  CLL125  CLL1464  CLL1323  CLL181  CLL115  CLL122  CLL373  CLL1292  CLL64  CLL60  CLL339  CLL824  CLL887 CLL5  CLL12  CLL16  CLL20  CLL23  CLL25  CLL26  CLL27  CLL30  CLL48  CLL56  CLL58  CLL63  CLL82  CLL83  CLL134  CLL138  CLL139  CLL141  CLL145  CLL157  CLL166  CLL176  CLL184  CLL188  CLL277  CLL282  CLL290  CLL294  CLL296  CLL306  CLL318  CLL356  CLL371  CLL386  CLL435  CLL442  CLL723  CLL749  CLL761  CLL853  CLL875  CLL1059  CLL1076  CLL1101  CLL1169  CLL1183  CLL1192  CLL1297  CLL1344  CLL1371  CLL1384  CLL1431  CLL148  CLL802  CLL1432  CLL10  CLL1191  CLL1430  CLL177  CLL358  CLL577  CLL1267  CLL15  CLL6  CLL677  CLL278  CLL128  CLL44  CLL753  CLL1078  CLL9  CLL199  CLL618  CLL880  CLL1446  CLL174  CLL372  CLL661  CLL1319  CLL29  CLL283  CLL1164  CLL367  CLL795  CLL804  CLL1424  CLL192  CLL467  CLL1179  CLL564  CLL745  CLL343  CLL393  CLL39  CLL832  CLL519  CLL1525  CLL783  CLL1360  CLL137  CLL90  CLL342  CLL33  CLL4  CLL654  CLL1353  CLL288  CLL1426  CLL1425  CLL1427  CLL84  CLL1276  CLL151  CLL308  CLL776  CLL244  CLL477  CLL32  CLL862  CLL1056  CLL628  CLL1423  CLL179  CLL473  CLL594  CLL523  CLL684  CLL1477  CLL125  CLL1464  CLL1323  CLL181  CLL115  CLL122  CLL373  CLL1292  CLL64  CLL60  CLL339  CLL824  CLL887  CLL1339 Sig.  6,7 framework (NNMF) extracted different mutational processes. Signature 9 (nc-AID) was extracted also among U-CLL in contrast with their known pathogenesis (b). This is a typical example of inter-sample bleeding and it can be solved either running a fitting approach after the initial NNMF analysis using only the catalog of signatures extracted by NNMF (c), or analyzing M-CLL and U-CLLs in two different and independent runs (d, e). Using the 30 COSMIC signatures as reference, the first approach is usually the most appropriate in order to estimate the real contribution of each single mutational process. In fact, the NNMF extracted signatures may be over or under split, therefore preventing a precise estimation of their contribution. For example, in this analysis, Signature 1 and 5 were extracted as one single process and only by running a fitting approach we were able to differentiate these two processes (c). Sig. = signature. In b and c, red patient labels are used for U-CLL, green for M-CLL, and blue for unknown cases

[C>A]A T[C>A]C T[C>A]G T[C>A]T A[C>G]A A[C>G]C A[C>G]G A[C>G]T C[C>G]A C[C>G]C C[C>G]G C[C>G]T G[C>G]A G[C>G]C G[C>G]G G[C>G]T T[C>G]A T[C>G]C T[C>G]G T[C>G]T A[C>T]A A[C>T]C A[C>T]G A[C>T]T C[C>T]A C[C>T]C C[C>T]G C[C>T]T G[C>T]A G[C>T]C G[C>T]G G[C>T]T T[C>T]A T[C>T]C T[C>T]G T[C>T]T A[T>A]A A[T>A]C A[T>A]G A[T>A]T C[T>A]A C[T>A]C C[T>A]G C[T>A]T G[T>A]A G[T>A]C G[T>A]G G[T>A]T T[T>A]A T[T>A]C T[T>A]G T[T>A]T A[T>C]A A[T>C]C A[T>C]G A[T>C]T C[T>C]A C[T>C]C C[T>C]G C[T>C]T G[T>C]A G[T>C]C G[T>C]G G[T>C]T T[T>C]A T[T>C]C T[T>C]G T[T>C]T A[T>G]A A[T>G]C A[T>G]G A[T>G]T C[T>G]A C[T>G]C C[T>G]G C[T>G]T G[T>G]A G[T>G]C G[T>G]G G[T>G]T T[T>G]A T[T>G]C T[T>G]G T[T>G
In general, our preferred approach to investigate mutational signatures in hematological malignancies follows three different steps: (1) signature discovery with a de novo extraction process; (2) assignment of extracted signatures to a reference catalog (i.e., COSMIC) and possibly identification of novel ones; (3) a fitting approach including only the subset of COSMIC signatures identified from the extraction process (Fig. 8). This multi-step approach allows the identification of known and novel signatures and their correct quantification, avoiding artefactual calls related to bleeding and overfitting. Based on a similar approach, two novel robust and stringent tools have recently been developed allowing the identification of >30 new mutational signatures and the redefinition of the previous 30-COSMIC signatures, creating a catalog to be used as reference for future studies 31 . These improved knowledge banks and bio-informatic tools will further refine our ability to investigate mutational signatures in hematological malignancies. However, we are convinced that prior knowledge of cancer biology and genomics will always be indispensable for a correct data interpretation.

Methods
Sample selection and processing of genomic data. In this study, we analyzed the single-nucleotide variant (SNV) catalog from four WGS cohorts: 143 CLLs (EGAS00000000092) 52,53 , 30 MMs (EGAD00001003309) 3,9 , 50 AMLs (phs000178. v1.p1) 59 , and two unpublished t-AML (EGAD00001005028). These last two cases were sequenced after written informed consent was obtained at the Wellcome Sanger Institute using the X10 Illumina platform. FASTQ files were aligned to the reference genome using BWAmem, and deduplicated aligned BAM files were analyzed using the following tools: ASCAT for copy number changes, BRASS for structural variations (large inversions and deletions, translocations, internal tandem duplication), Caveman and Pindel for Single-Nucleotide Variants (SNVs) and small insertion-deletions 20,66-68 , respectively. The characterization of the main clinical and genomic features of MM and CLL series is summarized in Supplementary Table 1 and Supplementary Data 2, respectively. Kataegis was defined as a cluster of 6 or more consecutive mutations with an average intermutation distance of less than or equal to 1 Kb 20 .
The study involved the use of human samples, which were collected after written informed consent was obtained (Wellcome Trust Sanger Institute protocol number 15/046 for the myeloma samples, Fondazione IRCCS Istituto Nazionale dei Tumori code 127/16 for the t-AML samples).
Mutational signature workflow. Mutational signatures were investigated using three published and available algorithms: the Alexandrov et al. 6 NNMF framework, deconstructSigs 24 and mutationalPatterns 37 R packages. The full mutatio-nalPatterns analysis was written in R and the code is provided in Supplementary Sofware Files 1-3 for MM, CLL, and AML respectively. Each of the above methods produces a matrix decomposition C ≈ SE, where C is the catalog matrix, with mutation types as rows and samples as columns, S is the signature matrix, with mutation types as rows and signatures as columns, and E is the exposure matrix, with signatures as rows and samples as columns ( Supplementary Fig. 1). The reconstruction error indicates how similar the mutational profiles of samples in C are to those in the product SE, and can be computed using different metrics, such as cosine similarity, Kullback-Leibler divergence (KLD) or RMSE.
Each of the signatures extracted with either mutationalPatterns or the method from Alexandrov et al. 6,7,37 were assigned to one or a combination of two COSMIC signatures. To do so, cosine similarities between the extracted signatures and each COSMIC signature, or a linear combination of two COSMIC signatures (using non-negative least squares R package NNLS), were computed. These results are available in Supplementary Data 1.
HRDetect in multiple myeloma. Analysis of homologous recombination deficiency (HRD) from BRCA1/BRCA2 deficiency as a possible source of genomic instability was performed using the recently published HRDetect algorithm 18 . The structural variant and indel catalog in MM were generated using BRASS and Pindel, respectively 20,67 .
Single-nucleotide variants on IGH. The mutation cancer cell fraction for c-AID SNVs was estimated using the Dirichlet process for both CLLs and MMs 4,9 . Considering the well-known complexity and low-quality mappping of IGH region, we ran three additional SNV callers (mutect2 69 , caveman 66 , and muse 70 ) to reduce the rate of false positives and we combined the results with the published catalog of SNVs generated with Sidron 52 . Seventy-nine percent of the previously published mutations on IGH was confirmed by at least one additional caller ( Supplementary  Fig. 12). Furthermore 512 additional SNVs were called by at least two out of the three new callers. Only mutations called by at least 2 out of 4 callers were included in the final analysis.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Code availability
All R codes used to generate signature data using mutationalPatterns in the paper are provided as Supplementary Software Files 1-3. All codes have been generated using R software v. 3.4.2.

Data availability
The sequencing data pertaining to MM are available from the European Genomephenome archive (EGA) database under the accession code EGAD00001003309. The sequencing data pertaining to CLL are available from EGA under the accession code EGAS00000000092. The published and unpublished AML sequencing data are available from dbGAP under the accession code phs000178 and from EGA dbGAP under the accession code EGAD00001005028, respectively. The breast cancer WGSs are available from the EGA under the accession code EGAS00001001178 20 . All the other data supporting the findings of this study are available within the article and its supplementary information files and from the corresponding author upon reasonable request. A reporting summary for this article is available as a Supplementary Information file.