Evaluation in heterogeneous cohorts

We now move to synthetic datasets with empirical heterogeneous signature weights. Here, we use absolute signature contributions (i.e., the number of mutations attributed to a signature) in WGS-sequenced samples from various cancers as provided by the COSMIC website (see Methods and Supplementary Fig. ⁶). For further evaluation, we choose eight types of cancer with mutually distinct signature profiles (Supplementary Fig. ⁷): Hepatocellular Carcinoma (Liver-HCC), Stomach Adenocarcinoma (Stomach-AdenoCA), Head and Neck Squamous Cell Carcinoma (Head-SCC), Colorectal Carcinoma (ColoRect-AdenoCA), Lung Adenocarcinoma (Lung-AdenoCA), Cutaneous Melanoma (Skin-Melanoma), Non-Hodgkin Lymphoma (Lymph-BNHL), and Glioblastoma (CNS-GBM). The first four cancer types all have SBS5 as the main contributing signature but substantially differ in the subsequent signatures. The remaining four cancer types have different strongly contributing signatures: SBS4, SBS7a, SBS5 and SBS40, and SBS40, respectively. The relative signature weights in individual samples were then used to construct synthetic datasets with a given number of mutations, allowing us to assess the performance of the fitting tools in realistic settings. Compared with the previous scenario with single-signature samples, there are now two main differences. First, nearly all samples have more than one active signature (the highest number of active signatures in one sample is eleven). Second, signature contributions differ from one sample to another; the average cosine distance between signature contributions in different samples ranges from 0.19 for Liver-HCC to 0.50 for ColoRect-AdenoCA. This scenario is thus referred to as heterogeneous cohorts.

We evaluate the fitting tools on heterogeneous cohorts with different numbers of mutations (100, 2,000, and 50,000) that cover the common range of mutational burden in WES and WGS analyses. Heterogeneous cohorts are more difficult to fit than the previously single-signature cohorts. For 2,000 mutations, for example, the lowest mean fitting error is 0.055 (achieved by SigProfilerAssignment) whereas the fitting error of mmsig for the same number of mutations is below 0.01 for all signatures. This is a direct consequence of heterogeneous signature weights: Even when the number of mutations is as high as 10,000, a signature with a relative weight of 2% contributes only 200 mutations and, as we have seen, the fitting errors are high for such a small number of mutations. Overall, SigProfilerSingleSample has the lowest fitting error for 100 mutations per sample. SigProfilerAssignment becomes the best method, with mmsig close second, for 2,000 mutations per sample. SigProfilerAssignment and MuSiCal are the best methods by a large margin for 50,000 mutations per sample. The best-performing methods are similar when binary classification metrics (precision, sensitivity, and the F1 score) are considered (see Supplementary Fig. ⁸ for further evaluation metrics). Fitting methods based on nonnegative least squares (represented by MutationalPatterns in Fig. 2) suffer from low precision due to overfitting. Flat signatures such as SBS5 and SBS40 are commonly underrepresented by the fitting tools (Supplementary Fig. ⁹) and the tools disagree on them most (Supplementary Fig. ¹⁰).

Open in a separate window

Fig. 2

A comparison of signature fitting methods on heterogeneous cohorts.

Mean fitting error (a) and precision and sensitivity (b) for different numbers of mutations per sample (columns) for the evaluated fitting tools. The results are averaged over 50 cohorts with 100 samples for eight distinct cancer types (see Methods). Each tool used all 67 COSMICv3 signatures as the reference catalog. The best-performing tool in each panel is marked with a frame (in the bottom row, the best tool by the F1 score combining precision and sensitivity; the dashed contours correspond to F1=0.9, 0.8, …). Results are not shown for SigsPack, YAPSA, and all three variants of sigminer as they are close (fitting error correlation above 0.999) to the results of MutationalPatterns.

Tool performance differs greatly by cancer type (Supplementary Fig. ¹¹). For 2,000 mutations, three different tools achieve the lowest fitting error for individual cancer types (SigProfilerAssignment for four of them, mmsig for three, and SigProfilerSingleSample for one). For 50,000 mutations, MuSiCal and SigProfilerAssignment are the best methods for five and three cancer types, respectively (see Supplementary Fig. ¹² for the dependence on the number of mutations). The results are similar when other performance metrics are considered (Supplementary Fig. ¹³). The overall similarity between the fitting error values (Supplementary Fig. ¹⁴) shows two distinct clusters of tools: the first formed by signature.tools.lib and all tools directly based on non-negative least squares and the second including deconstructSigs, sigLASSO, and MuSiCal.

One of the best-performing tools, SigProfilerSingleSample, reports a sample reconstruction similarity score that is often used to remove the samples whose reconstruction score is low. Our results show that this score is strongly influenced by the number of mutations; its absolute value is thus a poor indicator of the fitting accuracy (Supplementary Fig. ¹⁵). The same is the case of SigProfilerAssignment (Supplementary Fig. ¹⁶).

Choosing the reference catalog

The chosen reference catalog of signatures can significantly impact the performance of fitting algorithms^4,24. When the newer COSMICv3.2 catalog of mutational signatures is used as a reference instead of COSMICv3, the fitting error increases for most methods (Supplementary Fig. ¹⁷) due to an increased number of signatures (from 67 to 78) which makes the methods more prone to overfitting. However, SigProfilerSingleSample, SigProfilerAssignment, mmsig, MuSiCal, and sigLASSO are robust to increasing the number of reference signatures when the number of mutations per sample is high. On the other hand, when the reference signatures are constrained to the signatures that have been previously observed for a given cancer type and to artifact signatures⁴, the fitting error decreases substantially for most methods (Supplementary Fig. ¹⁸). Methods that are robust to adding reference signatures benefit less from removing irrelevant signatures, which in turn diminishes differences between methods (Supplementary Fig. ¹⁹).

To better illustrate the effect of reducing the number of reference signatures, we use three selected methods (SigsPack which together with MutationalPatterns and sigminer is the most sensitive to the reference catalog, SigProfilerSingleSample, and sigLASSO) and classify their output to true positive signatures, false positive signatures that are relevant to a given cancer type and false positive signatures that are irrelevant to a given cancer type (Fig. 3a). Using the whole COSMICv3 as a reference (top row), a simple method (SigsPack) starts with 17% of relevant false positives and 56% of irrelevant false positives. For comparison, these numbers are only 21% and 27% for SigProfilerSingleSample and 8% and 23% for sigLASSO (which, furthermore, leaves 50% of the mutations unassigned). When only relevant signatures are used as reference (bottom row), SigsPack improves much more than the two other methods. This further demonstrates how simple methods are particularly sensitive to the reference catalog and overfitting. Another interesting observation is that while, regardless of the reference catalog, SigsPack and sigLASSO tend to zero false positives as the number of mutations increases, this is not the case for SigProfilerSingleSample (as shown in Supplementary Fig. ¹¹, SigProfilerSingleSample performs poorly for Stomach-AdenoCA used in Fig. 3a). SigProfilerSingleSample evidently has a powerful algorithm to infer the active signatures from few mutations, but it does not converge to true signature weights when the mutations are many.

Open in a separate window

Fig. 3

Practical consequences of increasing the mutational burden.

a As the number of mutations per sample increases, weights assigned to false positives decrease (relevant false positives: excess weights assigned to signatures active in the cohort; irrelevant false positives: weights assigned to signatures not active in the cohort). We used simulated input data with 100 samples and Stomach-AdenoCA signature weights. The reference signature catalog is COSMICv3 (top row) and the 18 signatures active in the Stomach-AdenoCA samples (bottom row). b Performance of four fitting methods in identifying systematic differences in mutational weights between two groups of samples (see Methods for details). The success rate is the fraction of 250 synthetic cohorts with artificially introduced differences in SBS40 weights between even- and odd-numbered samples where the estimated SBS40 weights differ significantly (Wilcoxon rank-sum test, p-value below 0.05). The error bars show the 95% confidence interval (Wilson score interval).

Nevertheless, relying on a pre-determined list of relevant signatures is problematic for a number of reasons. First, the lists of signatures active in specific cancers are likely to change over time. Four of the eight considered cancer types have cohorts with <100 patients, so adding more WGS-sequenced tissues to the catalog is likely to significantly expand the list of signatures that are active in them. Second, most tools have difficulty recognizing that the provided signature catalog is insufficient even when the number of mutations in a sample is very large (Supplementary Fig. ²⁰; we return to this observation below). The estimated signature activities can therefore change in the future when better analytical tools become available. Third, when the list of relevant signatures is obtained from the COSMIC website, we rely on results obtained with one specific tool and this tool can be biased. We thus employ a different approach, which is based on fitting signatures using the whole COSMICv3 catalog (step 1) and then constraining the reference catalog to the signatures that sufficiently occur in the obtained results (step 2, see Methods for details).

The two-step fitting process can substantially improve the performance of some methods in some cases (Supplementary Fig. ²¹). At the same time, it tends to be detrimental to well-performing methods when the number of mutations per sample is high. It is then better to refrain from an ad hoc procedure for choosing a subset of reference signatures and instead rely on statistical selection mechanisms built into the methods themselves. A different multi-step process to select reference signatures was proposed in ref. ³. Signatures are first extracted de novo, each extracted signature is then assigned to one or two known reference signatures (see Methods for details), and thus-identified reference signatures are then used for fitting. This approach can also improve the fitting performance (Supplementary Fig. ²²). For samples with 100 mutations, the best results in terms of the fitting error and F1 score are obtained using SigProfilerSingleSample and SigProfilerAssignment, respectively, combined with the former two-step process to trim the reference signatures. For samples with 2,000 or 50,000 mutations, it is best to use the complete COSMIC catalog as a reference.

Signature fitting for downstream analysis

We have focused so far on estimates of signature weights and their errors with respect to a well-defined ground truth. In many cases, however, the estimates are only important as input for further downstream analyses assessing the correlations of signature weights with some clinicopathological parameters. We present here a simplified example of such an analysis by creating synthetic cohorts with CNS-GBM signature weights where statistically significant differences in the weights of signature SBS40 between even- and odd-numbered samples are introduced (see Methods for details). Estimation errors, together with the actual magnitude of the effect and the cohort size, are crucial to the ability to detect a significant difference in the signature weights. We contrast four fitting tools: simple MutationalPatterns, SigProfilerSingleSample and SigProfilerAssignment which perform well in Fig. 2, and widely used signature.tools.lib. When the number of mutations is large (10,000 in Fig. 3b), all tools are sufficiently precise to identify a statistically significant difference between the two groups of samples in nearly all cohorts. By contrast, when the number of mutations is smaller, the choice of the fitting tool matters. At 1000 mutations per sample, SigProfilerSingleSample is successfull in >90% of cohorts, while the other tools perform worse, in particular when the cohort size is small. At 100 mutations per sample, SigProfilerSingleSample still has some statistical power to detect a difference in SBS40 activities between the two groups of samples whereas the other tools fail regardless of the cohort size. When the signature of interest is easier to fit than SBS40 used in Fig. 3b, the differences between fitting tools become smaller (Supplementary Fig. ²³).

Synthetic mutational catalogs

Absolute contributions of signatures to sequenced tissues from various cancer types were downloaded from the Catalogue Of Somatic Mutations In Cancer (COSMIC, https://cancer.sanger.ac.uk/signatures/sbs/, catalog version 3.2) on November 14, 2021 for all 46 non-artifact SBS signatures that have tissue contribution data available. These data included only samples with reconstruction accuracy above 0.90 and for each sample, all signatures that contribute at least 10 mutations. Using the provided unique sample identifiers, we compiled the relative signature contributions in individual WGS-sequenced samples that are stratified by the cancer type. We chose eight different cancer types for further analysis: Hepatocellular carcinoma (Liver-HCC, n=422), stomach adenocarcinoma (Stomach-AdenoCA, n=170), head and neck squamous cell carcinoma (Head-SCC, n=57), colorectal carcinoma (ColoRect-AdenoCA, n=72), lung adenocarcinoma (Lung-AdenoCA, n=62), cutaneous melanoma (Skin-Melanoma, n=280), non-Hodgkin Lymphoma (Lymph-BNHL, n=117), and glioblastoma (CNS-GBM, n=63). All signatures in the tissue contribution data are from the COSMICv3 catalog.

We created a cohort of 100 samples with m mutations for each cancer type as follows. We first chose 100 samples (with repetition) from all samples for which empirical signature weights were available. Denoting the relative weights of signature s for sample i (sample composition) as w_si and the relative weight of context c for signature s (signature profile) as ω_cs, the relative weight of context c for sample i is obtained as a weighted sum over all signatures,

We then generated a synthetic mutational catalog with these weights by distributing m mutations among the 96 contexts, while the probability that a mutation is assigned to context c in sample i is W_ci. This is mathematically equivalent to first choosing signature s for sample i with probability w_si and then choosing context c with probability ω_cs. As a result, the number of mutations in contexts are multinomially distributed with mean mW_ci for context c and sampe i. The mean number of mutations contributed to sample i by signature s is mw_si. When this number is smaller than 10, the signature cannot be correctly identified by the evaluated tools as we use the standard procedure of setting the weights of signatures that contribute less than ten mutations to zero. To not bias the evaluation, we: (1) For sample i and m mutations per sample, remove all signatures that have mw_si<10 and (2) normalize the weights of the remaining signatures to one.

The approach described above allows us to reproduce empirical signature weights in previously analyzed samples without resorting to assumptions such as a log-normal distribution of the number of mutations due to a given signature^14,18 or adding additional zeros to the Poisson distribution to reproduce signatures that are not active in a sample²⁰. The code for generating synthetic mutational catalogs and evaluating the signature fitting tools, SigFitTest, is available at https://github.com/8medom/SigFitTest.

For Fig. 3b, we created synthetic mutational catalogs with signature activities modeled after real CNS-GBM samples where systematic differences in the activity of signature SBS40 have been introduced. After generating sample weights as described in the previous paragraph, the relative weights of SBS40 were multiplied by 1.3 and divided by 1.3 for even- and odd-numbered samples, respectively. The relative weights of all signatures in each sample were then normalized to one. Signature SBS40 is active in nearly all CNS-GBM samples (Supplementary Fig. ⁷) and its median relative weight is close to 50%. By the described process, the median weights in the two groups of samples become 42% and 55%, respectively. In another analysis (Supplementary Fig. ²³), we introduced differences in the weights of SBS1 for CNS-GBM samples by multiplying and dividing it by 1.2 in the two respective groups of samples.

To generate samples with signatures that are absent in the reference catalog COSMICv3, we used the 12 signatures that have been added in COSMICv3.3.1 (SBS10c, SBS10d, SBS86, SBS87, SBS88, SBS89, SBS90, SBS91, SBS92, SBS93, SBS94, SBS95). For a sample with the total weight of out-of-reference signatures x, we first multiply the weights of all COSMICv3 signatures by 1 − x. We then choose two signatures from the list above at random and assign them weights 0.7x and 0.3x (for Fig. 4d) or 0.5x each (Fig. 4e where x is fixed to 0.2, i.e., out-of-reference signatures contribute 20% of all mutations in each sample).

Tools for signature fitting

The task of fitting known mutational signatures to a mutational catalog consists of finding the combination of signatures from a reference set that “best” matches the given catalog. Denoting the matrix with reference signatures as $\mathsf{R}$, where R_cs is the relative weight of context c for signature s, and the given normalized mutational catalog as $\mathsf{m}$, where m_ci is the fraction of mutations in context c in sample i, this amounts to solving

with respect to $\mathsf{w}$, where w_si is the relative weight of signature s in sample i. To allow for a unique solution, vectors representing weights of different signatures must be linearly independent. The number of reference signatures thus cannot be higher than the number of mutational contexts (96 contexts for the common SBS signatures). Equation (2) is thus an over-determined system of linear equations. Several fitting tools are therefore based on minimizing the difference between m and $\mathsf{R}\mathit{w}$ through non-negative least squares (as the signature weights cannot be negative) or quadratic programming. We evaluated several tools that belong to this class: MutationalPatterns v3.14.0²⁶, YAPSA v1.30.0²⁹, SigsPack v1.18.0³⁰, and sigminer v2.3.1³¹ which has three separate methods based on quadratic programming, non-linear least squares, and simulated annealing. We find that all these tools produce similar results.

Other tools use various iterative processes by which the provided set of reference signatures is gradually reduced by removing the signatures, for example, whose inclusion does not considerably improve the match between the observed and reconstructed mutational catalogs (or, opposite, signatures are gradually added as long as the reconstruction accuracy sufficiently improves). We evaluated deconstructSigs v1.9.0³², SigProfilerSingleSample v0.0.0.27¹⁴, SigProfilerAssignment v0.1.7³³, mmsig v0.0.0.9000²⁴, and signature.tools.lib v2.2.0¹⁹, that all belong to this category. The newest tool, SigProfilerAssignment, combines backward and forward iterative adjustment of the reference catalog and these steps are repeated until convergence.

Finally, sigLASSO v1.1 combines the data likelihood in a generative model with L1 regularization and prior knowledge³⁴. To allow for a fair comparison with other tools, we did not use prior knowledge when evaluating sigLASSO. In our experience, this tool sometimes fails to halt but starting it again with the same input data resolves the issue. A similar Bayesian framework is used by sigfit v2.2.0³⁵. MuSiCal v1.0.0 is a recent tool which uses likelihood-based sparse NNLS to fit signatures³⁶.

While most tools produce signature weights that sum to one when normalized by the number of mutations, signature.tools.lib and sigLASSO explicitly allow for unassigned mutations. We used standard parameter settings for all tools. The authors of sigfit recommend setting all signatures whose lower bounds of the estimated 95% confidence intervals are below 0.01 to zero. The authors of deconstructSigs recommend setting all signatures whose relative weights are below 0.06 to zero. All results shown for sigfit and deconstructSigs follow these recommendations. Similarly to the COSMIC database, activities of signatures that contribute <10 mutations are set to zero. YAPSA makes it possible to use signature-specific cut-offs. However, the tool has only one set of precomputed cut-offs whose derivation is not clearly documented and these cut-offs work poorly for some cancer types. We thus do not use signature-specific cut-offs for YAPSA.

Evaluation metrics

Fitting error quantifies the difference between the true relative signature weights (which are used to generate the input mutational catalogs) and the estimated relative signature weights (for tools that estimate absolute signature weights, those are first normalized by the number of mutations in each sample). Denoting the true and estimated relative weights of signature s in sample i as w_si and ${\tilde{w}}_{si}$, respectively, the absolute error is summed over all signatures,

and then averaged over all samples. The division by two is introduced for normalization purposes. The lowest fitting error, 0, is achieved when the estimated signature weights are exact for all signatures and all samples. The highest fitting error, 1, is achieved when the estimated signature weights sum to one for each sample but they are all false positives. For example, when a sample has 40% contribution of SBS1 and 60% of SBS4 but a tool estimates 20% contribution of SBS2 and 80% contribution of SBS13, the fitting error is 1. We further quantify the agreement between the true and estimated signature weights by computing their Pearson correlation. For each sample where at least three signatures have positive either true or estimated weights, we compute the Pearson correlation coefficient between these two vectors whilst excluding all signatures that are zero for both of them. The obtained values are averaged over all samples in the cohort.

Further common evaluation metrics are based on classifying the estimated signatures as true positives (when both ${\tilde{w}}_{si}$ and w_si are positive), false positives (when ${\tilde{w}}_{si}>0$ and w_si=0), true negatives (when ${\tilde{w}}_{si}=0$ and w_si=0), and false negatives (when ${\tilde{w}}_{si}=0$ and w_si>0)^18,20,33. False positive weight quantifies how much weight is assigned to the signatures that are not active in the corresponding samples. The value for sample i,

is averaged over all samples. The lower the value, the better. Precision quantifies the reliability of the identified active signatures. It is computed by averaging the number of true positives to all positives over all samples in the cohort. Sensitivity quantifies the ability to identify all active signatures. It is computed by averaging the number of true positives to the number of active signatures (i.e., the sum of true positives and the false negatives) over all samples. Higher sensitivity can be commonly achieved at the cost of lower precision³⁷. This is addressed by the F1 score which is the harmonic mean of precision and recall (we compute the F1 score for each sample and then average it over all samples). To achieve a good F1 score, high precision and sensitivity are necessary.

Selecting reference signatures

It has been argued that the removal of irrelevant reference signatures can improve the fitting results⁵. To assess this hypothesis, we test a two-step process where: (1) We fit the samples using all COSMICv3 signatures as a reference and (2) keep only the signatures whose relative weight exceeds threshold w₀ for at least 5 samples in our cohorts with 100 samples. We use thresholds w₀=0.1 (m=100), w₀=0.03 (m=2000) and w₀=0.01 (m=50,000) which correspond to absolute signature contributions 10, 60, and 500, respectively. This pruning of reference signatures is often beneficial (Supplementary Fig. ²¹). However, for a high number of mutations and well-performing methods (mmsig, MuSiCal, sigLASSO, SigProfilerAssignment, SigProfilerSingleSample), this two-step process increases the produced fitting errors (m=50,000 in Supplementary Fig. ²¹). When lower thresholds are used to keep signatures, this can be avoided but there is nevertheless no improvement for m=50,000 and the improvements for m=100 vanish (results not shown).

We also assess the approach based on the extraction of signatures de novo as suggested in ref. ³. For signature extraction, we use SignatureAnalyzed v0.0.7 (https://github.com/getzlab/SignatureAnalyzer) with L1 prior on both W and H and a Poisson objective function for optimization. As recommended in ref. ³, each extracted signature is then compared with individual signatures from COSMICv3 as well as linear combinations of two signatures from COSMICv3. The signature (or a pair of signatures) that yields the smallest cosine distance is then added to the list of trimmed reference signatures for the given input data. When the number of mutations is small (m=100), less than five signatures are selected, on average. When the number of mutations is large (m=50,000), 90% of the active signatures are selected, on average. This approach can improve the fitting results (Supplementary Fig. ²²), in particular when the number of mutations is small. However, a direct comparison of the results obtained by the two described approaches shows that the former process based solely on signature fitting produces the lowest fitting error/highest F1 score for m=100 when combined with SigProfilerSingleSample and SigProfilerAssignment, respectively. For more mutations, it is best to use SigProfilerAssignment, mmsig, or MuSical and use all COSMICv3 signatures as a reference.

The recent tool MuSiCal³⁶ has a so-called full pipeline where de novo signature discovery is followed by matching the extracted signatures to a given catalog and refitting the data using the matched catalog signatures. The matching procedure is principled and involves automatic optimization of model parameters. At the same time, MuSiCal’s full pipeline is prohibitively computationally demanding for large-scale evaluation in many different settings as we do here.

We finally note that any approach to select a subset of reference signatures based on the input mutational catalogs is affected by the cohort size. For a small cohort, de novo signature extraction is likely to produce inferior results which will affect the subsequent refitting. All results shown here apply to cohorts of 100 samples.

Real mutational catalogs

To complement the synthetic data, we used real mutational catalogs produced by the International Cancer Genome Consortium (ICGC). Their PCAWG datasets¹⁴ are available from the ICGC upon request. From 2,780 WGS PCAWG samples, we kept the 146 samples that had at least 50,000 single base substitutions. Selecting the samples with high mutational burden is motivated by Fig. 2 where small fitting errors are achieved for 50,000 mutations. Estimates of signature activities in these samples by SigProfilerAssignment, sigLASSO, and MuSiCal (that perform best in Fig. 2 for 50,000 mutations) were used to choose the samples where the three tools disagree the least. To quantify the disagreement, we used the total absolute difference in estimated signature weights between two methods, averaged over all three pairs of methods. To prevent choosing samples that are overly similar to each other, we required that the total absolute difference between the mutational profiles must be at least 0.5. The input mutational profiles as well as the reconstructed profiles corresponding to the estimated signature weights agree well for the four chosen samples (Supplementary Fig. ²⁶). The ground truth signature estimates were obtained by averaging the positive signature weights over the three fitting tools; signatures that have been found active by less than two tools are set to zero. Thus-obtained signature weights were then normalized to one. Active signatures in sample 1 are SBS2 (50%), SBS13 (49%), and SBS1 (1%). Active signatures in sample 2 are SBS7a (88%) and SBS7b (12%). Active signatures in sample 3 are SBS5 (37%), SBS2 (36%), SBS13 (26%), and SBS1 (1%). Active signatures in sample 4 are SBS7a (58%), SBS7b (30%), SBS7c (7%), and SBS7d (5%). To study correlations of signature estimtes with clinical parameters, we used the clinical data of the OV-AU project contained in the ICGC Data Portal data Release 28 (2019-11-26) and the mismatch repair status of PCAWG samples (proficient/deficitent) provided by ref. ³⁶ in Supplementary Table ⁴.

We thank Peter M. Degen, Yitzhak Zimmer, Daniel M. Aebersold, and Kjong-Van Lehmann for their valuable support.