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.

sca_metadata = scaMetadata()

print(sca_metadata)

             Cancer_Type        Center NCBI_Build Sequence_Source
   gbm_tcga          GBM broad.mit.edu         37             WXS
   hnsc_tcga        HNSC broad.mit.edu         37         Capture
   kirc_tcga        KIRC broad.mit.edu         37         Capture
   luad_tcga        LUAD broad.mit.edu         37             WXS
   lusc_tcga        LUSC broad.mit.edu         37             WXS
   ov_tcga            OV broad.mit.edu         37             WXS
   skcm_tcga        SKCM broad.mit.edu         37         Capture
   thca_tcga        THCA broad.mit.edu         37             WXS
             Sequencing_Phase      Sequencer Number_Samples
   gbm_tcga           Phase_I Illumina GAIIx            291
   hnsc_tcga          Phase_I Illumina GAIIx            319
   kirc_tcga          Phase_I Illumina GAIIx            297
   luad_tcga          Phase_I Illumina GAIIx            538
   lusc_tcga          Phase_I Illumina GAIIx            178
   ov_tcga            Phase_I Illumina GAIIx            142
   skcm_tcga          Phase_I Illumina GAIIx            266
   thca_tcga          Phase_I Illumina GAIIx            406
             Number_Patients                           Cancer_Name
   gbm_tcga              291               Glioblastoma multiforme
   hnsc_tcga             319 Head and Neck squamous cell carcinoma
   kirc_tcga             293                    Kidney Chromophobe
   luad_tcga             519                   Lung adenocarcinoma
   lusc_tcga             178          Lung squamous cell carcinoma
   ov_tcga               142     Ovarian serous cystadenocarcinoma
   skcm_tcga             264               Skin Cutaneous Melanoma
   thca_tcga             403                    Thyroid carcinoma

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 [5]. Since the genomic positions are given in the NCBI notation and the references used later are in UCSC notation, the functions ucsc and ncbi are used to easily switch between the two notations. In this example, all mutational calls of a study will be pooled together, in order to find signatures related to a specific cancer type.

sca_vr = scaSNVRanges()

head(sca_vr, 3)

   VRanges object with 3 ranges and 1 metadata column:
         seqnames           ranges strand         ref              alt
                        
     gbm     chr1 [887446, 887446]      +           G                A
     gbm     chr1 [909247, 909247]      +           C                T
     gbm     chr1 [978952, 978952]      +           C                T
             totalDepth       refDepth       altDepth   sampleNames
            
     gbm                                   TCGA-06-5858
     gbm                                   TCGA-32-1977
     gbm                                   TCGA-06-0237
         softFilterMatrix |    study
                  | 
     gbm                  |      gbm
     gbm                  |      gbm
     gbm                  |      gbm
     -------
     seqinfo: 22 sequences from an unspecified genome
     hardFilters: NULL

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

sort(table(sca_vr$study), decreasing = TRUE)

   
     luad   skcm   hnsc   lusc   kirc    gbm   thca     ov 
   208724 200589  67125  61485  24158  19938   6716   5872

In a first step, the sequence motif for each variant is extracted based on the reference sequence. Here, the BSgenomes object for the human hg19 reference is used. However, all objects with a defined getSeq method can serve as the reference, e.g. an indexed FASTA file. Additionally, we transform all motifs to have a pyrimidine base (C or T) as a reference base [14].

sca_motifs = mutationContext(sca_vr, BSgenome.Hsapiens.UCSC.hg19, unify = TRUE)

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.

sca_mm = motifMatrix(sca_motifs, group = "study", normalize = TRUE)

head(round(sca_mm, 4))

             gbm   hnsc   kirc   luad   lusc     ov   skcm   thca
   CA A.A 0.0083 0.0098 0.0126 0.0200 0.0165 0.0126 0.0014 0.0077
   CA A.C 0.0093 0.0082 0.0121 0.0217 0.0156 0.0192 0.0009 0.0068
   CA A.G 0.0026 0.0061 0.0046 0.0144 0.0121 0.0060 0.0004 0.0048
   CA A.T 0.0057 0.0051 0.0070 0.0134 0.0100 0.0092 0.0007 0.0067
   CA C.A 0.0075 0.0143 0.0215 0.0414 0.0390 0.0128 0.0060 0.0112
   CA C.C 0.0075 0.0111 0.0138 0.0415 0.0275 0.0143 0.0018 0.0063

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.

plotMutationSpectrum(sca_motifs, "study")

Mutation spectrum over 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.

n_sigs = 5

sigs_nmf = identifySignatures(sca_mm, n_sigs, nmfDecomposition)

sigs_pca = identifySignatures(sca_mm, n_sigs, pcaDecomposition)

The results contains the decomposed matrices stored in a list and can be accessed using standard R accessor functions.

sigs_nmf

   MutationalSignatures:
     Samples (8): gbm, hnsc, ..., skcm, thca
     Signatures (5): S1, S2, S3, S4, S5
     Motifs (96): CA A.A, CA A.C, ..., TG T.G, TG T.T

sigs_pca

   MutationalSignatures:
     Samples (8): gbm, hnsc, ..., skcm, thca
     Signatures (5): S1, S2, S3, S4, S5
     Motifs (96): CA A.A, CA A.C, ..., TG T.G, TG T.T

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 also be modified 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.

plotSignatureMap(sigs_nmf) + ggtitle("Somatic Signatures: NMF - Heatmap")

Composition of somatic signatures estimated with the NMF, represented as a heatmap.

plotSignatures(sigs_nmf) + ggtitle("Somatic Signatures: NMF - Barchart")

Composition of somatic signatures estimated with the NMF, represented as a barchart.

plotObservedSpectrum(sigs_nmf)

plotFittedSpectrum(sigs_nmf)

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.

plotSampleMap(sigs_nmf)

Occurrence of signatures estimated with the NMF, represented as a heatmap.

plotSamples(sigs_nmf)

Occurrence of signatures estimated with the NMF, represented as a barchart.

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.

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 [3] [8] and could be adapted to the setting of somatic signatures with existing implementations (the sva and leapp packages). While this correction is not performed here, we exemplify the usage by taking the different sequencing technologies of the studies into account.

library(sva)

df = as(sca_metadata, "data.frame") ## sample x covariable
pheno = data.frame(s = unlist(df[ ,"Sequence_Source"]), c = unlist(df[ ,"Cancer_Type"]))
rownames(pheno) = gsub("(.*)_.*", "\\1", rownames(pheno))
mod = model.matrix(~ s + c, data = pheno)
mod0 = model.matrix(~ c, data = pheno)

sv = sva(sca_mm, mod, mod0, method = "irw")

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. 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 chromosome 1, based on a sample of 100'000 locations. Analogously, the k-mer occurrence across the human exome can be obtained easily.

k = 3
n = 1e5
chrs = "chr1"

chr1_ranges = as(seqinfo(BSgenome.Hsapiens.UCSC.hg19), "GRanges")
chr1_ranges = keepSeqlevels(chr1_ranges, chrs)

k3_chr1 = kmerFrequency(BSgenome.Hsapiens.UCSC.hg19, n, k, chr1_ranges)

k3_chr1

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.

head(sca_mm)

data(kmers)
norms = k3wg / k3we
head(norms)

sca_norm = normalizeMotifs(sca_mm, norms)

head(sca_norm)

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

citation("SomaticSignatures")

   
     Julian Gehring, Bernd Fischer, Michael Lawrence, Wolfgang
     Huber (2014): SomaticSignatures: Inferring Mutational
     Signatures from Single Nucleotide Variants. bioRxiv
     preprint, http://dx.doi.org/10.1101/010686
   
   A BibTeX entry for LaTeX users is
   
     @Article{,
       title = {SomaticSignatures: Inferring Mutational Signatures from Single Nucleotide Variants.},
       author = {Julian Gehring and Bernd Fischer and Michael Lawrence and Wolfgang Huber},
       year = {2014},
       journal = {bioRxiv},
       doi = {10.1101/010686},
       url = {http://dx.doi.org/10.1101/010686},
     }

We welcome emails with 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 emails and follow-up questions should take place over the Bioconductor mailing list, which serves as a repository of information.
• The subject line should contain SomaticSignatures and a few words describing the problem. First search the Bioconductor mailing list, for past threads which might have answered your question.
• 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 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:

source("http://bioconductor.org/biocLite.R")
biocLite("SomaticSignatures")

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.

library(VariantAnnotation)

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:

vr = VRanges(
    seqnames = "chr1",
    ranges = IRanges(start = 1000, width = 1),
    ref = "A",
    alt = "C")

vr

   VRanges object with 1 range and 0 metadata columns:
         seqnames       ranges strand         ref              alt
                    
     [1]     chr1 [1000, 1000]      +           A                C
             totalDepth       refDepth       altDepth   sampleNames
            
     [1]                                           
         softFilterMatrix
                 
     [1]                 
     -------
     seqinfo: 1 sequence from an unspecified genome; no seqlengths
     hardFilters: NULL