Inferring Somatic Signatures from Single Nucleotide Variant Calls

The SomaticCancerAlterations package provides the somatic SNV calls for eight WES studies, each investigating a different cancer type. The metadata summarizes the biological and experimental settings of each study.

The starting point of the analysis is a VRanges object which describes the somatic variants in terms of their genomic locations as well as reference and alternative alleles. For more details about this class and how to construct it, please see the documentation of the VariantAnnotation package [12]. In this example, all mutational calls of a study will be pooled together, in order to find signatures related to a specific cancer type.

To get a first impression of the data, we count the number of somatic variants per study.

In a first step, the sequence motif for each variant is extracted based on the genomic sequence. Here, the BSgenomes object of the human hg19 reference is used for all samples. However, personalized genomes or other sources for sequences, for example an indexed FASTA file, can be used naturally. Additionally, we transform all motifs to have a pyrimidine base (C or T) as a reference base [2]. The resulting VRanges object then contains the new columns context and alteration which specify the sequence content and the base substitution.

To continue with the estimation of the somatic signatures, the matrix \(M\) of the form {motifs × studies} is constructed. The normalize argument specifies that frequencies rather than the actual counts are returned.

The observed occurrence of the motifs, also termed somatic spectrum, can be visualized across studies, which gives a first impression of the data. The distribution of the motifs clearly varies between the studies.

The somatic signatures can be estimated with each of the statistical methods implemented in the package. Here, we will use the NMF and PCA, and compare the results. Prior to the estimation, the number \(r\) of signatures to obtain has to be fixed; in this example, the data will be decomposed into 5 signatures.

The individual matrices can be further inspected through the accessors signatures, samples, observed and fitted.

Up to now, we have performed the decomposition based on a known number \(r\) of signatures. In many settings, prior biological knowledge or complementing experiments may allow to determine \(r\) independently. If this is not the case, we can try to infer suitable values for \(r\) from the data.

Using the assessNumberSignatures function, we can compute the residuals sum of squares (RSS) and the explained variance between the observed \(M\) and fitted \(WH\) mutational spectrum for different numbers of signatures. These measures are generally applicable to all kinds of decomposition methods, and can aid in choosing a likely number of signatures. The usage and arguments are analogous to the identifySignatures function.

The obtained statistics can further be visualized with the plotNumberSignatures. For each tested number of signatures, black crosses indicate the results of individual runs, while the red dot represents the average over all respective runs. Please note that having multiple runs is only relevant for randomly seeded decomposition methods, as the NMF in our example.

\(r\) can then be chosen such that increasing the number of signatures does not yield a significantly better approximation of the data, i.e. that the RSS and the explained variance do not change sufficiently for more complex models. The first inflection point of the RSS curve has also been proposed as a measure for the number of features in this context [9]. Judging from both statistics for our dataset, a total of 5 signatures seems to explain the characteristics of the observed mutational spectrum well. In practice, a combination of a statistical assessment paired with biological knowledge about the nature of the data will allow for the most reliable interpretation of the results.

To explore the results for the TCGA data set, we will use the plotting functions. All figures are generated with the ggplot2 package, and thus, their properties and appearances can directly be modified, even at a later stage.

Focusing on the results of the NMF first, the five somatic signatures (named S1 to S5) can be visualized either as a heatmap or as a barchart.

Each signature represents different properties of the somatic spectrum observed in the data. While signature S1 is mainly characterized by selective C>T alterations, others as S4 and S5 show a broad distribution across the motifs.

In addition, the contribution of the signatures in each study can be represented with the same sets of plots. Signature S1 and S3 are strongly represented in the GBM and SKCM study, respectively. Other signatures show a weaker association with a single cancer type.

In the same way as before, the results of the PCA can be visualized. In contrast to the NMF, the signatures also contain negative values, indicating the depletion of a somatic motif.

Comparing the results of the two methods, we can see similar characteristics between the sets of signatures, for example S1 of the NMF and S2 of the PCA.

Since the observed mutational spectrum is defined by the data alone, it is identical for both all decomposition methods.

An alternative approach to interpreting the mutational spectrum by decomposition is clustering. With the clusterSpectrum function, the clustering is computed, by grouping either by the sample or motif dimension of the spectrum. By default, the Euclidean distance is used; other distance measures, as for example cosine similarity, are implemented is the proxy package and can be passed as an optional argument.

When investigating somatic signatures between samples from different studies, corrections for technical confounding factors should be considered. In our use case of the TCGA WES studies, this is of minor influence due to similar sequencing technology and variant calling methods across the studies. Approaches for the identification of so termed batch effects have been proposed [10,14] and existing implementations can be used in identifying confounding variables as well as correcting for them. The best strategy in addressing technical effects depends strongly on the experimental design; we recommend reading the respective literature and software documentation for finding an optimal solution in complex settings.

From the metadata of the TCGA studies, we have noticed that two different sequencing approaches have been employed, constituting a potential technical batch effect. The ComBat function of the sva package allows us to adjust for this covariate, which yields a mutational spectrum corrected for contributions related to sequencing technology. We can then continue with the identification of somatic signatures as we have seen before.

If comparisons are performed across samples or studies with different capture targets, for example by comparing whole exome with whole genome sequencing, further corrections for the frequency of sequence motifs can be taken into account [11]. The kmerFrequency function provides the basis for calculating the occurrence of k-mers over a set of ranges of a reference sequence.

As an example, we compute the frequency of 3-mers for the human toplevel chromosomes, based on a sample of 10'000 locations.

Analogously, the k-mer occurrence across a set of enriched regions, such as in exome or targeted sequencing, can be obtained easily. The following outlines how to apply the approach to the human exome.

With the normalizeMotifs function, the frequency of motifs can be adjusted. Here, we will transform our results of the TCGA WES studies to have the same motif distribution as of a whole-genome analysis. The kmers dataset contains the estimated frequency of 3-mers across the human genome and exome.

When we determine the sequence context for each alteration, we typically use one of the reference BSgenome packages in Bioconductor. But we are not restricted to those, and derive the somatic motifs from different types of sequence sources, for example 2bit and FASTA files. More precisely, the mutationContext function will work on any object for which a getSeq method is defined. You can get the full list available on your system, the results may vary depending on which packages you have loaded.

This allows us to perform our analysis also on non-standard organisms and genomes, for which a BSgenome package is not available, for example the 1000genomes human reference sequence. Or we can generate genomic references for specific populations, by updating the standard genomes with a set of known variants; see the documentation of the BSgenome package and the injectSNPs function in particular for this.

Taking further, we can base our analysis on the personalized genomic sequence for each individual, in case it is available. If we imagined that we had a set of somatic variant calls as VCF files and the personalized genomic sequence as FASTA files for two individuals A and B at hand, here a simple outline on how our analysis could work.

If you use the SomaticSignatures package in your work, please cite it:

We welcome questions or suggestions about our software, and want to ensure that we eliminate issues if and when they appear. We have a few requests to optimize the process:

• All questions and follow-ups should take place over the Bioconductor support site, which serves as a repository of information. First search the site for past threads which might have answered your question.
• The subject line should contain SomaticSignatures and a few words describing the problem.
• If you have a question about the behavior of a function, read the sections of the manual page for this function by typing a question mark and the function name, e.g. ?mutationContext. Additionally, read through the vignette to understand the interplay between different functions of the package. We spend a lot of time documenting individual functions and the exact steps that the software is performing.
• Include all of your R code and its output related to the question you are asking.
• Include complete warning or error messages, and conclude your message with the full output of sessionInfo().

Before you want to install the SomaticSignatures package, please ensure that you have the latest version of R and Bioconductor installed. For details on this, please have a look at the help packages for R and Bioconductor. Then you can open R and run the following commands which will install the latest release version of SomaticSignatures:

Over time, the packages may also receive updates with bug fixes. These installed packages can be updated with:

A central object in the workflow of SomaticSignatures is the VRanges class which is part of the VariantAnnotation package. It builds upon the commonly used GRanges class of the GenomicRanges package. Essentially, each row represents a variant in terms of its genomic location as well as its reference and alternative allele.

There are multiple ways of converting its own variant calls into a VRanges object. One can for example import them from a VCF file with the readVcf function or employ the readMutect function for importing variant calls from the MuTect caller directly. Further, one can also construct it from any other format in the form of: