Abstract

Background

Throughout their lives, cells are exposed to many different influences that can compromise the integrity of their DNA by introducing changes to the genome [1]. These somatic mutations are randomly introduced into the genome by various biochemical processes. These processes have different affinities for local genomic sequences, so that they leave their mark in the form of a specific mutation pattern on the genome of the cancer cell [2]. These patterns are termed (somatic) mutational signatures. Alexandrov et al. [3] have extracted and characterized a first set of mutational signatures and more signatures are being reported continuously.

The catalogue of somatic mutations in cancer, short COSMIC [4], hosts various sets of consensus signatures that have been found during analysis of the aggregation of multiple datasets from distinct types of human cancer [4, 5]. Some of these mutational signatures have been linked to environmental factors, like tobacco smoking or UV-light and the constituting DNA repair-mechanisms [3, 6], others have been associated with intrinsic processes such as defective DNA mismatch repair [7, 8]. The detection of these signatures in a tumor sample can thus yield helpful insights about the cancer’s aetiology [3] for diagnosis, prevention [9] and therapy [10].

Different frameworks have been proposed to mathematically decipher whether provided reference signatures are present in sequencing data from a single patient and how much they each contributed to its mutational load [11, 12]. Following Alexandrov et al. [3], the mutations are defined as the 6 single nucleotide variants (C >A:G >T, C >G:G >C, C >T:G >A, T >A:A >T, T >C:A >G & T >G:A >C) flanked by one nucleotide on each side. The tri-nucleotide formed by the mutated nucleotide and its two neighbours is called the context of the mutation. Mathematically, the mutational profile derived from sample data (mutational catalogue), can be expressed as vector m of K=96 somatic mutation frequencies that have been observed in the cancer sample. Mutational signatures are described by a matrix P which elements P_kn reflect the frequency with which the mutational process corresponding to the n^th signature causes the k^th mutational feature [3]. As the exposure of a mutational profile to a signature represents the signature’s contribution to the mutational load of the former, exposures to a set of processes can be inferred from a mutational catalogue by minimizing the difference between the observed & reconstructed catalogues:

The reconstructed catalogue Pe is the product between the mutational profiles matrix P and the exposure vector e, when individual exposures e_n are restricted to non-negative values. This formalism assumes that different processes have additive contributions to the mutational load.

We describe SigsPack, a Bioconductor package to estimate exposures to processes described by a known mutational signature matrix, for example from COSMIC. SigsPack also provides estimates of exposure stability, using bootstrapping. Its performance is benchmarked against synthetic & real data, using multiple tumor samples collected from the same patient. The effect of mutation context frequency is discussed, as well as the stability of individual COSMIC signatures and the loss of accuracy suffered by small mutational catalogues.

Implementation

Synthetic datasets

We used our package to construct two datasets consisting of simulated mutational catalogues using the function create_mut_catalogues. Samples in set1 are drawn from a distribution formed with equal contributions of signatures 7, 13, 21, 24 and 28 (so with each a weight of 0.2), likewise, the mutational catalogues of set2 have been samples from a distribution of signatures 3, 5, 8, 16 and 25 each having a weight of 0.2. The number of mutations in the catalogues were set to 1000, except when assessing the catalogue’s size effect on reconstructed exposures accuracy (see Fig. 3).

Open in a separate window

Fig. 3

Comparison of the mean prediction error with regard to the amount of mutations in the profile. Simulated samples have been created from set1 & set2 with different number of mutations. All have been bootstrapped 1000 times. The plot shows the mean prediction error, i.e. the mean for all SSE between the original exposure matrix and the predicted one, as a function of the mutational catalogue size (shown on logarithmic scale)

Results

Similarity between signatures

By construction, COSMIC signatures are non-orthogonal, in the sense that all the 96 mutation contexts are found in more than one signature. This non-orthogonality can affect the stability of the exposures. Figure 1 illustrates different aspects of the non-orthogonality between COSMIC signatures. Figure 1a displays the distance between signatures expressed as pairwise cosine similarity. Non-orthogonality effects can be assessed by computing the error between each signature profile and its reconstruction using all the 29 other signatures. Figure 1b shows that four signatures (5, 6, 19 and 26) have a cosine similarity higher than 0.95 with their reconstructed profile. Those signatures might be labelled as unstable, as the information contained in their profile is mostly contained in the remaining 29 others.

Open in a separate window

Fig. 1

Similarity between COSMIC signatures. a Pairwise cosine similarity of the 30 COSMIC signatures. Co-linear signatures have a similarity of 1, and orthogonal signatures a similarity of 0. Signatures 5 & 16 have a similarity higher than 0.9. b Computed cosine similarity of each signature with the profile that was constructed by trying to reconstruct each signature with the 29 others. The signatures are ordered by decreasing similarity to their reconstructed profile

Accuracy of exposure estimation

To investigate the possible effects of the non-orthogonality between COSMIC signatures, we have measured the accuracy of exposure reconstruction using synthetic datasets (see the “Implementation” section). Based on Fig. 1, we have selected 5 fairly distant signatures (set1, signatures 7, 13, 21, 24 & 27), and 5 relatively similar signatures (set2, signatures 3, 5, 8, 16 & 25), and generated mutational catalogues from each of them. Figure 2a shows that in certain cases, exposures can be faithfully recapitulated, and that exposures are relatively robust to small changes in the mutational catalogue. Inferred exposures for signatures absent from set1 are very small (Additional file 2): the 3rd quartile of the bootstrapped values for those signatures is above 1% only for signatures 1 & 17, and the bootrapped maximum values above 10% for signatures 4, 16, 17 & 29. As signature estimations cannot be negative by construction, the overall effect of estimation errors is to contribute to a small underestimation of the contribution of 4 out of 5 present signatures (7, 13, 21 & 27).

Open in a separate window

Fig. 2

Signature prediction on a simulated mutational catalogue. The catalogue is sampled from set1 signatures (7, 13, 21, 24 and 27), each occurring with a probability of 20% (a), and from set2 (3, 5, 8, 16 and 25), with the same 20% occurring probability (b). In both case, the mutation catalogue consists of 1000 mutations and was bootstrapped 1000 times. The signature contribution was predicted with quadratic programming, for each re-sample the distribution is shown. The stars mark the prediction for the original profile

The reconstruction of exposures in mutational catalogues drawn from set2 (Fig. 2b & Additional file 2) shows that any small amount of noise in the data leads to dramatic changes in the exposure predictions. Exposure to signature 5 is underestimated in more than 75% of the re-sampled catalogues, and half re-sampled catalogues underestimate signature 25 contribution by almost 50%. Signature 26 (which is absent from the generation protocol) is assigned a contribution higher than 3.6% in 25% of the re-sampled catalogues.

Required number of mutations

In order to quantify how the precision of the inferred exposures depends on the mutational load of a sample, simulated mutational catalogues based on set1 and set2 have been created with different amounts of mutations. Figure 3 displays the mean prediction error, i.e. the mean difference between the predicted signatures exposure and the actual one. In general, we observe that the mean prediction error is higher on samples with only a small mutational load, consistent with Rosenthal et al. [12]. We also notice that the prediction error is lower on set1, constructed with stable signatures, than it is on set2, constructed with unstable signatures, even for high mutational loads.

Bootstrapping can faithfully recover the biological variance

Next, we have investigated the stability of exposures inference in the presence of biological variability, using the 13 samples described in the “Implementation” section. As the experimental protocol was identical for all samples, the differences between inferred inferences within the same patient could then be attributed mostly to the heterogeneities in the clonal compositions of each sample. Figure 4 shows two examples of the agreement between exposure reconstruction from different samples. In that figure, the bootstrapping estimates have been computed from a single sample (the complete set of plots can be found in Additional file 7).

Open in a separate window

Fig. 4

Exposures estimations for multiple samples from the same patient. a Patient 99, variability intervals obtained by bootstrapping sample 4. b Patient 69, variability intervals obtained by bootstrapping sample 4. In both cases, 1000 re-sampling realisations of each catalogue have been used by bootstrapping

If we postulate that boostrapping can provide estimates to exposure variability, then the bootstrapping intervals computed from the different samples should largely overlap, and the overlap should contain the original exposures’ estimation values. To test this hypothesis, we have computed, for each patient and each signature, the intersection of the exposure range obtained by bootstrapping from all samples. We have then asked whether any exposure estimated for individual samples fell outside of this range. This happens only for signature 27 (patients 69 & 80) and 5 (patient 99). In two out of three cases, the exposure estimation is higher than 0 (0.2% & 6%) for one sample, while the range is limited to 0 for another sample. In the third case, (signature 27 in patient 80), the exposure is estimated to be 2.49% in sample 2, while the bootstrapping range computed from sample 3 is between 0 and 2.43%. We conclude from these examples that in most cases, the boostrapping variability estimates can provide a proxy for the exposures’ variability due to clonal heterogeneity of mutational catalogues.

Across all patients, the contribution from unstable signatures vary considerably from one sample to the next, and there is more than 10% difference between estimates from the 4 samples for signatures 5, 9 & 16 for patient 69, and for signatures 3 & 16 for patient 99 (the corresponding signatures for patient 80 are 3, 5, 8 & 16). As for synthetic data, bootstrapping provides large variation estimate for those signatures: for signatures 5 & 16, the extreme values reached by bootstrapping differ by more than 30% for 10 out of the 13 samples. These difference appear to be compensating: across the 13 samples, the correlation between signatures 3 & 16 exposures is −0.75, the most negative correlation between all signature pairs. This observation supports the hypothesis that the unstable mutations can be explained almost equally well by their contributions. The complete set of exposures is found in Additional file 4. Signatures 6 & 10, often observed in colorectal cancers, are not present at high levels in any of the 13 samples analysed here. However, these signatures are associated with microsatellite instability (MSI) and mutations in POLE and the DNA repair mechanism [17]. Neither MSI nor POLE mutations were present in any of the 3 patients.

An example of the sequence composition effect

COSMIC signatures have been normalized to the trinucleotide frequencies of the human reference genome version GRCh37 [18]. The regions on which the mutations catalogues have been observed directly affect the exposures estimation (see Additional file 1 for details). In particular, in presence of mutational catalogues generated from exome data, the scaling of the mutational catalogue to the genome tri-nucelotide frequencies in not equivalent to the scaling of mutational signatures to the exome tri-nucleotide frequencies. Figure 5 illustrates the sequence composition effect for the first sample of patient 69 (the full set of exposure estimations after both normalization are found in Additional file 6). Although the frequencies of enriched regions and of the whole genome are reasonably similar (correlation coefficient above 0.72, Additional file 5), there is a difference greater than 10% in the inferred exposures of signature 16. While this might be attributed to the unstability of signature 16, signature 1 also shows a difference of 9.4% in exposure estimation, due to normalization. Over all samples and all patients, signature 1 appears to be most sensitive to the scaling choice: the estimation differs by more than 10% in 11 of the 13 samples, even more often than unstable signatures (signature 16 estimation is affected in 4 samples). As signature 1 is associated with endogenous C to T mutations, its contribution is mainly affected by the tri-nucleotides containing a central C. These are more frequent (relatively) in the exome regions than in the whole genome (Wilcoxon test P value 4.7·10⁻⁴). This simple example shows that taking sequence composition into consideration by normalisation of mutational catalogues is required to ensure accurate exposure values.

Open in a separate window

Fig. 5

Exposures estimations for different normalization of the same sample. a Normalization on the genome, where the observed mutation frequencies were re-scaled. To perform the scaling, tri-nucleotide ratios where first obtained by dividing the tri-nucleotide frequencies observed in exome regions by their corresponding values in the whole genome. The mutational catalogue was then scaled by the tri-nucleotide ratio and converted to frequencies. b Normalization on the exome, where the COSMIC signatures were re-scaled. In this case, the COSMIC signatures were scaled by the inverse of the tri-nucleotide ratios, and converted to frequencies again. Together, a) & b) show that the scaling direction (mutational catalogue or signatures) lead to different exposures estimations. In both cases, 1000 realisations have been used by bootstrapping

Exposures from random catalogues

To gain understanding on possible causes of the unstability of estimation of some signatures, we have created 1000 mutational catalogues of 1000 events each, drawn with probability proportional to the underlying tri-nucleotide frequencies in the human genome. Exposures were computed from these mutational catalogues, and the sum of estimated exposures over the 1000 samples is shown in Fig. 6. These "null" catalogues represent the absence of mutational process signal, as the occurrence frequencies reflect the corresponding frequencies in the genome. Estimated exposures from signatures 3 & 9 are higher than 50% and 25% respectively for more than half of the catalogues. However, 5 or 16 are consistently absent of the exposure estimation from most "null" catalogues, unlike the unstability displayed in presence of signal. In this respect, the lack of stability in the exposure estimation cannot be solely attributed to the truly random component of the mutational catalogue, which frequency approximates the tri-nucleotides frequencies in the genome.

Open in a separate window

Fig. 6

Estimated exposures from 1000 random mutational catalogues. These catalogues have been randomly drawn according to the human genome tri-nucleotide frequencies. Each catalogue has 1000 mutations

Discussion

Using both synthetic data & multiple samples from the same donors, we have shown that stability of exposure estimates can be accurately represented by bootstrapping mutational catalogues. It should be noted that reliable exposure estimation can only be achieved for a relatively large mutational catalogue. When the number of available mutations falls below a few hundreds, the number of observations in each of the 96 mutation classes is not sufficient to ensure a stable exposure reconstruction. SigsPack will warn the user when the mutational catalogue size is below 125.

Even for larger mutational catalogues, some exposures displayed considerable variability upon re-sampling of the mutational catalogue. Inspection of relationships between signatures suggest that similarity between signatures, and the fact that some signatures can be approximated by others, may cause unwanted sensitivity to mutational catalogue details for some signatures’ exposures.

To overcome the problem of unstable signatures, it may be possible to follow for example Letouzé and co-workers [19] and use only a subset of signatures known to be involved in a specific cancer entity. Careful selection of signatures would presumably reduce considerably the redundancy shown in the COSMIC set, and decrease the exposures variability. The signature selection might also be guided by the particulars of the cohort under consideration. Also, for some mutational processes, the analysis of the presence of specific di-nucleotide mutations and/or indel in the patient’s somatic mutations might be more informative than the quantification of the associated signatures [10]. Such signatures might possibly be omitted, again leading to a decrease of exposures variability.

We have shown that accurate exposures estimation requires matching tri-nucleotide frequencies between regions on which mutational catalogues and signature matrix have been collected. When using the COSMIC signature matrix provided with SigsPack, exome data must be put on “genomic scale” prior to exposures estimation. Choice of reference tri-nucleotide frequencies and normalization should be carefully selected when creating new mutational signature matrix from a large cohort. Even though SigsPack doesn’t provide algorithms to generate such matrix, it allows for easy normalization of each of its component. We suggest that scaling observations to the whole genome’s tri-nucleotide frequencies should be encouraged, as it renders the results independent of the experimental particulars.

The analysis above has been carried out on the full set of COSMIC version 2 signatures. However, as the recently released version 3 contains more signatures, the signatures will remain non-orthogonal, and there might be cluster of signatures very similar to each other. These features could presumably lead to the same stability problems for exposure estimation as with COSMIC version 2 signatures.

Conclusions

Many computational methods aimed at mutational signatures discovery already exist [20]. The SigsPack package is aimed at estimating exposure to known mutational signatures, rather than the process of uncovering new ones. It relies on the Alexandrov additive mutational frequency model, rather than a position weight matrix [21] or a probablistic framework, such as EMu [22] or sigfit [23]. It builds on existing methods ([11–13]), and provides support for vcf input, exome regions normalization, exposure estimation and stability estimates for individual signatures. It also provides facilities to plot and examine estimated exposures, and a data generation module to benchmark user’s defined signature matrices.

Availability and requirements

• Project name SigsPack

• Project home page https://github.com/bihealth/SigsPack

• Operating system(s) any OS running R [24] & Bioconductor [25]

• Programming language R

• Other requirements Bioconductor

• License GPL-3

• Any restriction to use by non-academics None

Additional files

Acknowledgements

Abbreviations

Authors’ contributions

DB conceived and designed the project. FS performed the analysis & wrote the package. CM & EB contributed to the interpretation of results. TB & AB provided the samples. FS, DB, CM & EB wrote the manuscript. All authors have read and approved the manuscript.

Funding

Partially supported by Deutsche Forschungsgemeinschaft (SFB-TR36; T.B.), Deutsche Krebshilfe (111546; T.B.) and the Berlin Institute of Health (CRG-1; T.B.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The software package has been submitted to Bioconductor ([25]) and is available from github (https://github.com/bihealth/SigsPack). The patients mutational catalogues are included in Additional file 3. The tri-nucleotides frequencies for normalization are included in Additional file 5.

Ethics approval and consent to participate

Ethic approval was obtained by the local ethic committe of the Charité - Universitätsmedizin Berlin (E A1/265/14). All patients gave written informed consent for sample acquisition and analyses of patient data.

Consent for publication

All patients gave written informed consent for publication of patient data.

Competing interests

The authors declare that they have no competing interests.

Footnotes

Contributor Information

Franziska Schumann, Email: [email protected].

Eric Blanc, Email: [email protected].

Clemens Messerschmidt, Email: [email protected].

Thomas Blankenstein, Email: ed.nilreb-cdm@eknalbt.

Antonia Busse, Email: [email protected].

Dieter Beule, Email: [email protected].

References