A Sticky Multinomial Mixture Model of Strand-Coordinated Mutational Processes in Cancer.

The characterization of mutational processes in terms of their signatures of activity relies mostly on the assumption that mutations in a given cancer genome are independent of one another. Recently, it was discovered that certain segments of mutations, termed processive groups, occur on the same DNA strand and are generated by a single process or signature. Here we provide a first probabilistic model of mutational signatures that accounts for their observed stickiness and strand coordination. The model conditions on the observed strand for each mutation and allows the same signature to generate a run of mutations. It can both use known signatures or learn new ones. We show that this model provides a more accurate description of the properties of mutagenic processes than independent-mutation achieving substantially higher likelihood on held-out data. We apply this model to characterize the processivity of mutagenic processes across multiple types of cancer.

Introduction

Mutational processes are key factors in shaping cancer genomes (Alexandrov et al., 2013a, Alexandrov et al., 2013b, Helleday et al., 2014, Tubbs and Nussenzweig, 2017), and their characterization has important implications for understanding the disease and choosing targeted therapies (Davies et al., 2017a, Davies et al., 2017b, Polak et al., 2017, Gulhan et al., 2019). Multiple algebraic and statistical approaches have been suggested for the detection of mutational processes from somatic mutation data (Alexandrov et al., 2013a, Alexandrov et al., 2013b, Fischer et al., 2013, Shiraishi et al., 2015, Kim et al., 2016, Rosales et al., 2016). These methods, which focus on single-base substitutions (consult Figure 1), are based on learning the pattern of mutations of each potential process as well as its activity (aka exposure) in any given tumor in a way that will best explain the observed mutation data. State-of-the-art approaches for learning mutational signatures include non-negative matrix factorization (NMF) methods (Alexandrov et al., 2013a, Alexandrov et al., 2013b, Fischer et al., 2013, Kim et al., 2016, Rosales et al., 2016) that aim to explain the mutation counts as a sum over all signatures of the probability of a specific mutation to be generated by the respective signature times its exposure. Other approaches that borrow from the world of topic modeling (discovery of abstract topics in text documents) aim to provide a probabilistic model of the data so as to maximize the model's likelihood (Shiraishi et al., 2015, Funnell et al., 2018). However, most of these methods assume that mutations are independent of one another and cannot capture processes that create dependencies among them.

Figure 1

Definitions and Conventions

The figure shows normal DNA, mutated DNA, representation and characteristics of DNA mutations, and different types of stickiness used in StickySig model variants.

(A) The genome consists of the reference strand (the strand whose 5′-end is on the short arm of the chromosome), also known as the Watson strand or the plus strand, and the complementary strand, also known as the Crick strand or the minus strand.

(B) In the mutated DNA, changes in DNA base pair sequence are shown in red.

(C) Single base-pair substitution (SBS) can be represented as a transition in which one base pair in normal DNA is replaced by another in mutated DNA. The reference allele refers to the nucleotide that is found in the reference strand of normal genomic DNA. The pyrimidine strand is the strand (+or –) containing the pyrimidine base (C or T) in normal DNA. It is usually not known which base in a pair was the source of a mutation; thus, the convention is to annotate mutations from the pyrimidine base of the mutated base pair, leading to 6 substitution types (when context is not considered) or to 96 possible combinations of substitution types and neighboring bases.

(D) The StickySig model can use several types of stickiness opportunities: all mutations can be sticky, same strand stickiness, mutations having a pyrimidine base in the normal DNA on the same strand as the previous mutation; same allele stickiness, mutations having the same reference allele as the previous mutation; same substitution stickiness, mutations having exactly the same base-pair substitution as the previous mutation; same mutation stickiness, mutations having the same mutation features (96 mutation category and pyrimidine strand) as the previous mutation; and none, no stickiness allowed that leads to MMM model. Other types of stickiness can be also considered.

Recently, it was observed that some signatures operate in a strand-coordinated manner where consecutive mutations tend to mutate from the same reference allele and occur on the same strand (Nik-Zainal et al., 2014). Morganella et al. generalized these observations and found segments of such mutations (i.e., same reference allele and same strand) that they termed processive groups (Morganella et al., 2016). The length of a processive group, that is, the number of such consecutive mutations attributed to the same signature, is signature dependent. The significance and abundance of these processive groups suggested that certain mutational processes display stickiness and strand-coordination properties. In a previous work we have suggested a hidden Markov-based model for capturing sequential dependencies between close-by mutations (Wojtowicz et al., 2019). Here we follow a similar path and suggest novel probabilistic models for consecutive, albeit not necessarily close-by, mutations that occur on the same strand.

The biological reasons for this strand coordination are related, at least in part, to the asymmetric role that the two strands play in many cellular processes that operate on DNA. For example, the APOBEC C-to-U editing enzymes are a major source of mutations in many cancer types and are known to act on single-stranded DNA (Refsland and Harris, 2013). Many cellular processes, including replication and transcription, require strand separation leaving one or both strands exposed. Importantly, if the strands are separated, one of the strands is often more exposed than the other, leading to asymmetric strand coordination of APOBEC mutations. In particular, during DNA replication, the two DNA strands are processed differently. In this process, one of the strands (the lagging strand) is more exposed than the other strand (the leading strand). Owing to these differences, APOBEC asymmetry between these two strands is particularly strong (Haradhvala et al., 2016, Morganella et al., 2016, Tomkova et al., 2018, Seplyarskiy et al., 2016). In addition, leading and lagging strands are, among other differences, also processed by different polymerases, which might introduce different types of error in each strand leading to replication-related strand coordination. Transcription-coupled repair is another source of strand-specific mutagenesis. Another process leading to coordinated mutations and strand asymmetry is the formation DNA/RNA duplexes—the so-called R-loops. R-loops are thought to form co-transcriptionally when nascent messenger RNA hybridizes with the DNA template and thus can protect this strand from APOBEC activity and other types of mutations that act on single-stranded DNA. Indeed, multiple signatures have been found to have mutation strand bias in template versus non-template strands (Alexandrov et al., 2013a, Alexandrov et al., 2013b, Haradhvala et al., 2016, Morganella et al., 2016, Tomkova et al., 2018).

Our suggested probabilistic model, StickySig, accounts for the stickiness and strand coordination of mutational signatures. The model captures independent mutations as well as processive groups in one probabilistic framework. In cross-validation tests on multiple datasets, StickySig outperforms independent-mutation models or sticky models that do not account for the strand information. We apply our model to gain new insights about the stickiness and exposures of known signatures, as well as in a de novo setting to learn new signatures.

Results and Discussion

Mutation Data

For each dataset, we followed the standard approach introduced by Alexandrov et al., 2013a, Alexandrov et al., 2013b and classified mutations into M = 96 categories based on the 5′ flanking base, substitution type, and 3′ flanking base, following the convention that simple base-pair substitutions can be classified into six subtypes (Figures 1A–1C). (It is usually not known which base in a pair was the source of a mutation, thus the convention is to annotate substitutions from the pyrimidine base, i.e., G:C > A:T is written as C > T rather than G > A.) These mutations are assumed to be the result of the activity of K mutational processes, each of which is associated with a signature S=(e(1), …,e(M)) of probabilities to emit each of the mutation categories. Henceforth, we denote the mutation categories observed in a given tumor by o1, …,o. We assume that o was emitted by signature s (whose identity is hidden from us).

In addition to mutation categories, the mutation data include information on the reference allele (the nucleotide that is on the Watson strand) and the pyrimidine strand (the Watson or Crick strand containing the pyrimidine base in normal DNA). By convention, the Watson strand is the reference genome strand (the strand whose 5′-end is on the short arm of the chromosome) and the Crick strand is the complementary strand. See Figures 1A–1C.

Data Description

We analyzed breast cancer (BRCA), chronic lymphocytic leukemia (CLLE), and malignant lymphoma (MALY) mutation datasets from whole-genome sequences from the International Cancer Genome Consortium (ICGC) (more information is available in the Data and Code Availability supplement section). We chose to study BRCA, CLLE, and MALY because they are known to have active signatures (Signatures 2 and 13) that were previously shown to display strand coordination (Morganella et al., 2016). In addition, each of the corresponding datasets contained at least 100 samples.

A Comparative Evaluation

We evaluated our suggested models and compared them with previous approaches using the datasets outlined above (Table 1). Here MMM serves as a stand-in for state-of-the-art non-probabilistic mutation signature methods such as non-negative matrix factorization, as MMM is a probabilistic method that encodes the standard assumption that each mutation in a tumor is independent of all others. In the Supplemental Informationwe show that MMM is in fact equivalent to a statistical variant of NMF, which is widely used for mutational signature analysis (Fischer et al., 2013, Kim et al., 2016).

Table 1

Datasets Analyzed in This Study: Breast Cancer (BRCA), Malignant Lymphoma (MALY), and Chronic Lymphocytic Leukemia (CLLE)

Cancer Type	#Samples	#Mutations	COSMIC Signatures
BRCA	560	3,479,652	1, 2, 3, 5, 6, 8, 13, 17, 18, 20, 26, 30
MALY	100	1,220,526	1, 2, 5, 9, 13, 17
CLLE	100	270,870	1, 2, 5, 9, 13

We performed this comparison in two modes; the first was testing the models in a refitting mode when the signatures are known, and the second was testing the models in a de novo mode when the signatures are unknown. In the following comparisons, owing to running time considerations, we used a maximum of 100 iterations.

First, we focused on the refitting scenario, using the known COSMIC signatures. To this end, we use the leave-one-chromosome-out (LOCO) method. Specifically, we split samples by chromosomes and learned for each sample i the exposure vector π and the stickiness for the cosmic signatures α using all the chromosomes but one. We then report the log likelihood of the model on the left-out chromosome (summed across all samples and chromosomes). The results are summarized in Table 2 and clearly show the superiority of StickySig across the three cancer types analyzed. In each cancer type, StickySig has higher held-out likelihood than the independent mutation MMM, demonstrating that mutation signatures have stickiness that is shared across samples and that modeling this stickiness provides greater predictive power for held-out data. Furthermore, the difference between the variants of StickySig and MMM becomes much larger when StickySig is restricted to allow stickiness only between mutations with the same reference allele or the same base-pair substitution.

Table 2

Performance Evaluation of MMM and StickySig Variants in a Refitting Setting Using the Leave-One-Chromosome-Out (LOCO) Method

Model	LOCO Log likelihood
	BRCA	MALY	CLLE
MMM	−13743198	−5235042	−1178024
StickySig	−13739451	−5232119	−1177527
StickySig-same-strand	−13736711	−5233842	−1177905
StickySig-same-allele	−13696283	−5205381	−1173271
StickySig-same-substitution	−13549757	−5206208	−1173221
StickySig-same-mutation	−13683356	−5227289	−1176916

In bold are the best values for each dataset.

Next, we study the de novo scenario, where signatures are learned as part of the model training. For comparison purpose, we set the number of signatures to be the same as the number of active COSMIC signatures used above. To evaluate the models with respect to signature learning we used 10-fold sample cross-validation, where we learned e and α across all samples of the train data; then, for each sample X in the test set we fitted π and computed Pr[X|π,e,α] and summed across all samples. This tests the model's capability to produce meaningful signatures that can explain well a new given sample. The results are summarized in Table 3 and show again that stickiness adds power to the model. Here again the leading models are StickySig-same-allele and StickySig-same-substitution.

Table 3

Performance Evaluation of MMM and StickySig Variants in a De Novo Setting Using 10-Fold Sample Cross-Validation (CV)

Model	10-Fold CV Log likelihood
	BRCA	MALY	CLLE
MMM	−13713230	−5200582	−1167727
StickySig	−13708734	−5187470	−1167238
StickySig-same-strand	−13702999	−5202680	−1167628
StickySig-same-allele	−13231739	−4987710	−1121929
StickySig-same-substitution	−13135156	−5020921	−1128093
StickySig-same-mutation	−13597846	−5187859	−1165856

In bold are the best values for each dataset.

Finally, we wished to assess the signatures learned by the algorithm. We trained the two best performing variants of StickySig (StickySig-same-allele and StickySig-same-substitution) on each of the three datasets. For evaluation purpose, we matched each signature learned by each model to its most similar COSMIC signature known to be active in the corresponding cancer type (measured via cosine similarity). The results are summarized in Figure 2. Evidently, StickySig-same-allele yields signatures that are more similar to the COSMIC ones. Note that the definition of signatures in our model is different from the standard one since they are coupled with a stickiness value, which in this application was always one. This may partially explain the deviation (particularly in the same-substitution case) from the COSMIC signatures.

Figure 2

Signatures Learning

Performance of signature learning by StickySig-same-allele (A) and StickySig-same-substitution (B) on the three cancer datasets. For each case, depicted are the cosine similarities of the learned signatures to known COSMIC signatures, sorted from highest to lowest and computed by a maximum matching algorithm to prevent repetitions.

Strand-Coordinated StickySig Defines Processive Groups in Breast Cancers

Morganella et al. defined processive groups as sets of adjacent substitutions of the same mutational signature sharing the same reference allele (Morganella et al., 2016). Our model, StickySig-same-allele, allows us to compute maximum likelihood estimates that sequences of mutations are generated in processive groups; hence, we could apply it to characterize the processivity of the different signatures in breast cancer. In order to compare and contrast our findings with those of Morganella et al. (2016), we used the same statistical test for the significance of a processive group of a given length. We confirmed the association of processive groups with Signatures 2, 6, 13, 17, and 26 when using the same length threshold of more than 10 (Morganella et al., 2016). In addition, our strand-coordinated model revealed that processivity is also a feature of Signature 18 (Figure 3A). The number of processive groups of length more than 10 was particularly high for Signatures 2 and 13 (Figure 3B), which is consistent with previous studies showing that APOBEC-related signatures demonstrated strand-coordinated mutagenesis (Nik-Zainal et al., 2014, Morganella et al., 2016).

Figure 3

Stickiness in BRCA

(A) Relationship between processive group lengths (columns) and mutational signatures (rows) modeled by StickySig-same-allele. The size of each circle represents the number of groups (log10) observed for the specified group length and for each signature. The color of each circle corresponds to the p value of detecting a processive group of a given length in randomized data (-log10).

(B) The number of processive groups of length more than 10 for all signatures modeled by MMM (gray), StickySig (blue), and StickySig-same-allele (red).

(C) The total number of mutations can be sticky in StickySig (blue) and StickySig-same-allele (red).

(D) The total number of sticky mutations as modeled by StickySig (blue) and StickySig-same-allele (red).

(E) The number of sticky mutations for each signature as modeled by StickySig (blue) and StickySig-same-allele (red).

(F) Signature stickiness α as learned by StickySig (blue) and StickySig-same-allele (red). All bar plots show mean values with standard error of the mean (small black bars) from 10 random initializations of StickySig models.

Next, we tested whether using the strand-coordinated StickySig rather than the mutation-independent MMM or the regular StickySig was important for the accurate discovery of processive groups. On average, the StickySig-same-allele model uncovered 133.9 groups of length greater than 10 in 41.9 patients, whereas MMM and StickySig models captured only 38.5 in 11 patients and 38.3 in 11.9 patients, respectively (Figure 3B). This is consistent with the large number of sticky mutations found by StickySig-same-allele model, even though StickySig has more sticky opportunities (Figures 3C and 3D). All these differences underscore the higher sensitivity of the strand-coordinated model for detecting processive groups, which may in part explain the observed differences in likelihood between MMM and StickySig models on held-out data (Table 2).

Processive groups, as summarized in Figures 3A and 3B, capture statistically significant patterns in cancer genomes and are considered features of specific signatures. An alternative characterization signature could be provided by the model parameter α—the “stickiness” of a signature—which is learned by the strand-coordinated model. Thus, we analyzed how these two views of strand-coordinated mutagenesis relate to each other. We considered only signatures for which there is a sufficient number of sticky mutations to properly learn this parameter (Figure 3E). For comparison purposes, we also included stickiness values computed with the strand-oblivious StickySig. We found that, in the strand-coordinated variant, the most sticky signatures were 1, 2, 6, 13, 17, 18, 26, and 30 (Figure 3F). Signatures 2, 6, 13, 17, 18, and 26 are exactly the same signatures that were found to be associated with processive groups. Although Signature 30 did not make the length 10 cutoff for processive groups, its processive segments are also relatively long. Interestingly Signatures 3 and 5 were not found to be sticky despite the fact that their processive groups were also quite long. In contrast, there is some stickiness to Signature 1 while its processive groups are shorter than for 3, 5, and 30. The signature stickiness in the strand-oblivious StickySig model is minimal. The meaning of these intriguing findings is a subject for further investigation.

Limitations of the Study

Although in this work we focused on Watson/Crick strands, there is compelling evidence that other strand definitions may affect signatures and their activities. As we reviewed above, such categorization may be based on replication-based characteristics such as leading and lagging or transcription-based characteristics such as template and non-template. In that vein, a promising next step may be modeling multiple strand characteristics simultaneously, rather than considering them individually. For example, there is evidence in humans and other species that transcription and replication are co-oriented (Huvet et al., 2007, Srivatsan et al., 2010). Testing these different variants may reveal the role of strand characteristics in mutagenesis.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.