Introduction to MutationalPatterns

In MutationalPatterns: Comprehensive genome-wide analysis of mutational processes

knitr::opts_chunk$set(echo = TRUE)

options(width = 96)
library(ggplot2)
library(BiocStyle)

Introduction

Mutational processes leave characteristic footprints in genomic DNA. This package provides a comprehensive set of flexible functions that allows researchers to easily evaluate and visualize a multitude of mutational patterns in base substitution catalogues of e.g. healthy samples, tumour samples, or DNA-repair deficient cells. This is the second major version of the package. Many new functions have been added and functions from the previous version have been enhanced. The package covers a wide range of patterns including: mutational signatures, transcriptional and replicative strand bias, lesion segregation, genomic distribution and association with genomic features, which are collectively meaningful for studying the activity of mutational processes. The package works with single nucleotide variants (SNVs), insertions and deletions (Indels), double base substitutions (DBSs) and larger multi base substitutions (MBSs). The package provides functionalities for both extracting mutational signatures de novo and determining the contribution of previously identified mutational signatures on a single sample level. MutationalPatterns integrates with common R genomic analysis workflows and allows easy association with (publicly available) annotation data.

Background on the biological relevance of the different mutational patterns, a practical illustration of the package functionalities, comparison with similar tools and software packages and an elaborate discussion, are described in the new MutationalPatterns article, which is in currently being written. The old article can be found here.

This vignette shows some common ways in which the functions in this package can be used. It is however not exhaustive and won't show every argument of every function. You can view the documentation of a function by adding a ? in front of it. Like: ?plot_spectrum. The describes the functions and all their arguments in more detail. It also contains more examples of how the functions in this package can be used.

Installation

First you need to install the package from Bioconductor.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("MutationalPatterns")

You also need to load the package. This needs to be repeated every time you restart R.

library(MutationalPatterns)

Data

To perform the mutational pattern analyses, you need to load one or multiple VCF files with variant calls and the corresponding reference genome.

List reference genome

You can list available genomes using r Biocpkg("BSgenome"):

library(BSgenome)
head(available.genomes())

Make sure to install the reference genome that matches your VCFs. For the example data this is BSgenome.Hsapiens.UCSC.hg19.

Now you can load your reference genome:

ref_genome <- "BSgenome.Hsapiens.UCSC.hg19"
library(ref_genome, character.only = TRUE)

Load example data SNVs

We provided two example data sets with this package. One consists of a subset of somatic SNV catalogues of 9 normal human adult stem cells from 3 different tissues (Blokzijl et al. 2016), and the other contains somatic indels and DBSs from 3 healthy human hematopoietic stem cells (Osorio et al. 2018). The MBS data you will find in the latter dataset was manually included by us to demonstrate some features of this package.

This is how you can locate the VCF files of the example data from the first set.

These will be used for the SNV examples:

vcf_files <- list.files(system.file("extdata", package = "MutationalPatterns"),
  pattern = "sample.vcf", full.names = TRUE
)

You also need to define corresponding sample names for the VCF files:

sample_names <- c(
  "colon1", "colon2", "colon3",
  "intestine1", "intestine2", "intestine3",
  "liver1", "liver2", "liver3"
)

Now you can load the VCF files into a GRangesList:

grl <- read_vcfs_as_granges(vcf_files, sample_names, ref_genome)

Here we define relevant metadata on the samples, such as tissue type. This will be usefull later.

tissue <- c(rep("colon", 3), rep("intestine", 3), rep("liver", 3))

Load example data Indels, DBSs and MBSs

We will now locate the VCF files of the example data from the second set. These will be used for the Indels, DBS and MBS examples.

blood_vcf_fnames <- list.files(
  system.file("extdata", package = "MutationalPatterns"), 
  pattern = "blood.*vcf", full.names = TRUE)

Set their sample names.

blood_sample_names <- c("blood1", "blood2", "blood3")

Read in the data, without filtering for any mutation type using the type="all" argument. (By default only SNVs are loaded for backwards compatibility.)

blood_grl <- read_vcfs_as_granges(blood_vcf_fnames, blood_sample_names, 
                                  ref_genome, type = "all")

You can now retrieve different types of mutations from the GrangesList.

snv_grl <- get_mut_type(blood_grl, type = "snv")
indel_grl <- get_mut_type(blood_grl, type = "indel")
dbs_grl <- get_mut_type(blood_grl, type = "dbs")
mbs_grl <- get_mut_type(blood_grl, type = "mbs")

It's also possible to directly select for a specific mutation type when reading in the data. This can be a convenient shortcut, when you are only interested in a single type of mutation.

indel_grl <- read_vcfs_as_granges(blood_vcf_fnames, blood_sample_names, 
                                  ref_genome, type = "indel")

Mutation characteristics

SNVs

Base substitution types

You can retrieve base substitution types from the VCF GRanges object as "REF>ALT" using mutations_from_vcf:

muts <- mutations_from_vcf(grl[[1]])
head(muts, 12)

You can retrieve the base substitutions from the VCF GRanges object and convert them to the 6 types of base substitution types that are distinguished by convention: C>A, C>G, C>T, T>A, T>C, T>G. For example, when the reference allele is G and the alternative allele is T (G>T), mut_type returns the G:C>T:A mutation as a C>A mutation:

types <- mut_type(grl[[1]])
head(types, 12)

To retrieve the sequence context (one base upstream and one base downstream) of the base substitutions in the VCF object from the reference genome, you can use the mut_context function:

context <- mut_context(grl[[1]], ref_genome)
head(context, 12)

Withtype_context, you can retrieve the types and contexts for all positions in the VCF GRanges object. For the base substitutions that are converted to the conventional base substitution types, the reverse complement of the sequence context is returned.

type_context <- type_context(grl[[1]], ref_genome)
lapply(type_context, head, 12)

With mut_type_occurrences, you can count mutation type occurrences for all VCF objects in the GRangesList. For C>T mutations, a distinction is made between C>T at CpG sites and other sites, as deamination of methylated cytosine at CpG sites is a common mutational process. For this reason, the reference genome is needed for this functionality.

type_occurrences <- mut_type_occurrences(grl, ref_genome)
type_occurrences

Mutation spectrum

A mutation spectrum shows the relative contribution of each mutation type in the base substitution catalogs. The plot_spectrum function plots the mean relative contribution of each of the 6 base substitution types over all samples. Error bars indicate the 95% confidence interval over all samples. The total number of mutations is indicated.

p1 <- plot_spectrum(type_occurrences)

You can also plot the mutation spectrum with distinction between C>T at CpG sites and other sites:

p2 <- plot_spectrum(type_occurrences, CT = TRUE)

Other options include plotting the spectrum with the individual samples as points and removing the legend to save space:

p3 <- plot_spectrum(type_occurrences, CT = TRUE, 
                    indv_points = TRUE, legend = FALSE)

Here we use the r CRANpkg("gridExtra") package to combine the created plots, so you can see them next to each other.

library("gridExtra")
grid.arrange(p1, p2, p3, ncol = 3, widths = c(3, 3, 1.75))

It's also possible to create a facet per sample group, e.g. plot the spectrum for each tissue separately:

p4 <- plot_spectrum(type_occurrences, by = tissue, CT = TRUE, legend = TRUE)

Or you could use the standard deviation instead of a 95% confidence interval:

p5 <- plot_spectrum(type_occurrences, CT = TRUE, 
                    legend = TRUE, error_bars = "stdev")

grid.arrange(p4, p5, ncol = 2, widths = c(4, 2.3))

96 mutational profile

First you should make a 96 trinucleodide mutation count matrix. (In contrast to previous versions this also works for single samples.)

mut_mat <- mut_matrix(vcf_list = grl, ref_genome = ref_genome)
head(mut_mat)

Next, you can use this matrix to plot the 96 profile of samples. In this example we do this for 2 samples:

plot_96_profile(mut_mat[, c(1, 7)])

Larger contexts

It's also possible to look at larger mutational contexts. However, this is only usefull if you have a large number of mutations.

mut_mat_ext_context <- mut_matrix(grl, ref_genome, extension = 2)
head(mut_mat_ext_context)

The extension argument also works for the mut_context and type_context functions.

You can visualize this matrix with a heatmap.

plot_profile_heatmap(mut_mat_ext_context, by = tissue)

You can also visualize this with a riverplot.

plot_river(mut_mat_ext_context[,c(1,4)])

Indels

First you should get the COSMIC indel contexts. This is done with get_indel_context, which adds the columns muttype and muttype_sub to the GRangesList. The muttype column contains the main type of indel. The muttype_sub column shows the number of repeat units. For microhomology (mh) deletions the mh length is shown.

indel_grl <- get_indel_context(indel_grl, ref_genome)
head(indel_grl[[1]], n = 3)

Next count the number of indels per type. This results in a matrix that is similar to the mut_mat matrix.

indel_counts <- count_indel_contexts(indel_grl)
head(indel_counts)

Now you can plot the Indel spectra. The facets at the top show the indel types. First the C and T deletions. Then the C and T insertions. Next are the multi base deletions and insertions. Finally the deletions with microhomology are shown. The x-axis at the bottom shows the number of repeat units. For mh deletions the microhomology length is shown.

plot_indel_contexts(indel_counts, condensed = TRUE)

You can also choose to only plot the main contexts, without taking the number of repeat units or microhomology length into account.

plot_main_indel_contexts(indel_counts)

DBSs

First get the COSMIC DBS contexts. This is done by changing the REF and ALT columns of the GRangesList.

head(dbs_grl[[1]])
dbs_grl <- get_dbs_context(dbs_grl)
head(dbs_grl[[1]])

Next count the number of DBSs per type. This again results in a matrix that is similar to the mut_mat matrix.

dbs_counts <- count_dbs_contexts(dbs_grl)

Finally we can plot the DBS contexts. The facets at the top show the reference bases. The x-axis shows the alternative variants.

plot_dbs_contexts(dbs_counts, same_y = TRUE)

We can also choose to plot based on only the reference bases. Now the x-axis contains the reference bases.

plot_main_dbs_contexts(dbs_counts, same_y = TRUE)

MBSs

No COSMIC MBS contexts existed when this vignette was written. Therefore the length of the MBSs is used as its context. First we can count the MBSs. This again results in a matrix that is similar to the mut_mat matrix.

mbs_counts <- count_mbs_contexts(mbs_grl)

Next we can plot the contexts

plot_mbs_contexts(mbs_counts, same_y = TRUE)

Pooling samples

Sometimes you have very few mutations per sample. In these cases it might be usefull to combine multiple samples. This can be done with pool_mut_mat. This works on the matrixes of SNVS, Indels, DBSs and MBSs.

pooled_mut_mat <- pool_mut_mat(mut_mat, grouping = tissue)
head(pooled_mut_mat)

Mutational signatures

Mutational signatures are thought to represent mutational processes, and are characterized by a specific contribution of mutation types with a certain sequence context. Mutational signatures can be extracted de novo from your mutation count matrix, with non-negative matrix factorization (NMF). It's also possible to identify the exposure of your mutation count matrix to previously defined mutational signatures. This is often referred to as signature refitting. NMF is most usefull for large amounts of samples, while signature refitting can also be used on single samples. We will first discuss NMF and then signature refitting. Finally we will discuss analyzing the similarity between a mutational profile and signatures directly.

De novo mutational signature extraction using NMF

NMF

A critical parameter in NMF is the factorization rank, which is the number of mutational signatures you extract. You can determine the optimal factorization rank using the r CRANpkg("NMF") package [@Gaujoux2010]. As described in their paper:

``...a common way of deciding on the rank is to try different values, compute some quality measure of the results, and choose the best value according to this quality criteria. The most common approach is to choose the smallest rank for which cophenetic correlation coefficient starts decreasing. Another approach is to choose the rank for which the plot of the residual sum of squares (RSS) between the input matrix and its estimate shows an inflection point.''

In general, larger datasets allow you to use a higher rank. Here we will show NMF for SNVs. Performing NMF on other mutation types works the same way.

First add a small psuedocount to your mutation count matrix:

mut_mat <- mut_mat + 0.0001

Use the NMF package to generate an estimate rank plot. This can take a long time:

library("NMF")
estimate <- nmf(mut_mat, rank = 2:5, method = "brunet", 
                nrun = 10, seed = 123456, .opt = "v-p")

And plot it:

plot(estimate)

Extract mutational signatures from the mutation count matrix with extract_signatures. In this example 2 signatures are extracted, because a rank of 2 is used. (For larger datasets it is wise to perform more iterations by changing the nrun parameter to achieve stability and avoid bad local minima):

nmf_res <- extract_signatures(mut_mat, rank = 2, nrun = 10, single_core = TRUE)

NMF also works on other mutation types like Indels and DBS. You can even combine matrixes from different mutation types to, for example, extract combined Indel/DBS signatures.

combi_mat = rbind(indel_counts, dbs_counts)
nmf_res_combi <- extract_signatures(combi_mat, rank = 2, nrun = 10, single_core = TRUE)

Bayesian NMF

It's also possible to use variational bayes NMF. This could make it easier to determine, the correct rank. To do this you need to install the r Biocpkg("ccfindR") package. You can then determine the optimal number of signatures, which can again take a long time. Warnings will occur when you use ranks that are too high. (In this example we avoid these warnings by using nrun=1, combined with a set seed. In practice you shouldn't use a rank that's too high and you should also use a higher number for nrun.) With a larger dataset you could try higher ranks. The highest value in the plot is the mathematically optimal number of signatures. (A note of warning: The mathematically optimal number doesn't necessarily make biological sense.)

# BiocManager::install("ccfindR")
library("ccfindR")
sc <- scNMFSet(count = mut_mat)
set.seed(4)
estimate_bayes <- vb_factorize(sc, ranks = 1:3, nrun = 1, 
                               progress.bar = FALSE, verbose = 0)
plot(estimate_bayes)

Extracting the signatures is then done by:

nmf_res_bayes <- extract_signatures(mut_mat, rank = 2, nrun = 10, 
                                    nmf_type = "variational_bayes")

Changing the names of the extracted signatures

You can provide the extracted signatures with custom names:

colnames(nmf_res$signatures) <- c("Signature A", "Signature B")
rownames(nmf_res$contribution) <- c("Signature A", "Signature B")

It's possible that some of the signatures extracted by NMF are very similar to signatures that are already known. Therefore, it might be usefull to change the names of the NMF signatures to these already known signatures. This often makes it easier to interpret your results.

To do this you first need to read in some already existing signatures. Here we will use signatures from COSMIC (v3.1) [@Alexandrov2020]. (We will discuss how to use other signature matrixes later.)

signatures = get_known_signatures()

You can now change the names of the signatures extracted by NMF. In this example the name of a signature is changed if it has a cosine similarity of more than 0.85 with an existing COSMIC signature.

nmf_res <- rename_nmf_signatures(nmf_res, signatures, cutoff = 0.85)
colnames(nmf_res$signatures)

We now see that the signatures we extracted are very similar to COSMIC signatures SBS1 and SBS5. This helps with the interpretation because the aetiology of SBS1 is already known. This also tells us we didn't identify any completely novel processes. An extracted signature that is not similar to any previously defined signatures, is not proof of a "novel" signature. The extracted signature might be a combination of known signatures, that could not be split by NMF. This can happen when, for example, too few samples were used for the NMF.

Visualizing NMF results

You can plot the 96-profile of the signatures (When looking at SNVs):

plot_96_profile(nmf_res$signatures, condensed = TRUE)

You can visualize the contribution of the signatures in a barplot:

plot_contribution(nmf_res$contribution, nmf_res$signature,
  mode = "relative"
)

The relative contribution of each signature for each sample can also be plotted as a heatmap with plot_contribution_heatmap, which might be easier to interpret and compare than stacked barplots. The signatures and samples can be hierarchically clustered based on their euclidean distance. Clustering here is based on the similarity between the contributions. (Signatures with a similar contribution will thus be clustered together. The same applies for samples.)

Plot signature contribution as a heatmap with sample and signature clustering dendrograms:

plot_contribution_heatmap(nmf_res$contribution, 
                          cluster_samples = TRUE, 
                          cluster_sigs = TRUE)

It's also possible to provide your own signature and sample order. This can be a manual ordering, but in this example we use clustering. We can cluster the signatures based on their cosine similarity and then retrieve the order:

hclust_signatures <- cluster_signatures(nmf_res$signatures, method = "average")
signatures_order <- colnames(nmf_res$signatures)[hclust_signatures$order]
signatures_order

We can do the same thing for the samples:

hclust_samples <- cluster_signatures(mut_mat, method = "average")
samples_order <- colnames(mut_mat)[hclust_samples$order]
samples_order

Now we can use the signature and sample order in the contribution heatmap:

plot_contribution_heatmap(nmf_res$contribution,
  sig_order = signatures_order, sample_order = samples_order,
  cluster_sigs = FALSE, cluster_samples = FALSE
)

A reconstructed mutational profile has been made for each sample by the NMF, based on the signatures and their contribution. The better the NMF worked the more similar the reconstructed profile will be to the original profile.

We can compare the reconstructed mutational profile with the original mutational profile for a single sample like this:

plot_compare_profiles(mut_mat[, 1],
  nmf_res$reconstructed[, 1],
  profile_names = c("Original", "Reconstructed"),
  condensed = TRUE
)

This is the function for SNVs. For indels you would use plot_compare_indels, for DBSs, plot_compare_dbs and for MBSs plot_compare_mbs.

We can also plot the cosine similarity between the original and reconstructed matrix for all the samples. When a reconstructed profile has a cosine similarity of more than 0.95 with the original, the reconstructed profile is considered very good.

plot_original_vs_reconstructed(mut_mat, nmf_res$reconstructed, 
                               y_intercept = 0.95)

Signature refitting

Find mathematically optimal contribution of COSMIC signatures

Signature refitting quantifies the contribution of any set of signatures to the mutational profile of a sample. This is specifically useful for mutational signature analyses of small cohorts or individual samples, but also to relate own findings to known signatures and published findings. The fit_to_signatures function finds the optimal linear combination of mutational signatures that most closely reconstructs the mutation matrix by solving a non-negative least-squares constraints problem.

Fit mutation matrix to the COSMIC mutational signatures:

fit_res <- fit_to_signatures(mut_mat, signatures)

The fit_res object can be visualized similarly to the nmf_res object. The functions plot_contribution, plot_contribution_heatmap, plot_compare_profiles and plot_original_vs_reconstructed will all work. As an example we show the contribution of signatures as a barplot.

plot_contribution(fit_res$contribution,
  coord_flip = FALSE,
  mode = "absolute"
)

We also show the cosine similarity with the reconstructed profiles, as this gives a good idea of how well the signatures could explain the mutational profiles.

plot_original_vs_reconstructed(mut_mat, fit_res$reconstructed, 
                               y_intercept = 0.95)

Stricter refitting

In the previous plots, many signatures were used to explain the mutational profiles of the samples. It seems however unlikely that this many mutational processes were really active in these samples. This problem, known as overfitting, occurs because fit_to_signatures finds the optimal combination of signatures to reconstruct a profile. It will use a signature, even if it improves the fit very little.

Another issue with signature refitting is signature misattribution. Mutations will sometimes be attributed to different signatures in samples with a similar mutational profile. This can give the impression that samples are very different, when they actually aren't. This is often the result of "flat" signatures, which are harder to fit. Signatures that are similar to each other can also cause this misattribution issue.

One way to deal with overfitting and the misattribution of signatures is by selecting a limited number of signatures that will be used for the refitting. When you are analyzing a liver sample you could for example only use signatures that are known to occur in liver. This method is recommended by @Degasperi2020. Using prior knowledge like this will reduce overfitting, but can also introduce bias. You won't be able to identify signatures, if you removed them beforehand. Another downside of this method is that you need prior knowledge of which signatures could be present. We recommend using this method when possible.

Another way of dealing with overfitting is by starting with a standard refit and then removing signatures that have little effect on how well a mutational profile can be reconstructed. This works in an iterative fashion. In each iteration the signature with the lowest contribution is removed and refitting is repeated. Each time the cosine similarity between the original and reconstructed profile is calculated. You stop removing signatures when the difference between two iterations becomes bigger than a certain cutoff. This way only the signatures that are really necessary to explain a mutational profile will be used. This method is similar to a method used by @Alexandrov2020. In MutationalPatterns it can be used with fit_to_signatures_strict.

A downside of this method is that the cutoff you should use is somewhat subjective. Here we use a cutoff of 0.02. Decreasing this number will make the refitting less strict, while increasing it will make the refitting more strict.

strict_refit <- fit_to_signatures_strict(mut_mat, signatures, max_delta = 0.004)

This function returns a list containing a fit_res object and a list of figures, showing in what order signatures were removed during the refitting.

Here we show the figure for one sample. The x-axis shows the signature that was removed during that iteration. The red bar indicates that the difference in cosine similarity has become too large. The removal of signatures is stopped and SBS1 is kept for the final refit.

fig_list <- strict_refit$sim_decay_fig
fig_list[[1]]

The fit_res can be visualized in the same way as other fit_res objects.

fit_res_strict <- strict_refit$fit_res
plot_contribution(fit_res_strict$contribution,
  coord_flip = FALSE,
  mode = "absolute"
)

A third method that can reduce overfitting and the misattribution of signatures is to merge similar signatures. This works by merging signatures whose cosine similarity is higher than a certain cutoff value. These merged signatures can then be used for refitting. The benefit of this method is that you don't need prior knowledge. For most common use-cases, we don't recommend this method, because it is less conventional and can be harder to interpret. However, we provide it here to give you the possibility to use it if you need it. You can merge signatures like this:

merged_signatures <- merge_signatures(signatures, cos_sim_cutoff = 0.8)

The best refitting method will depend on your data and research question. A single method can be used, but it's also possible to combine several methods.

Bootstrapped refitting.

The stability of signature refitting can be suboptimal, because of the previously mentioned signature misattribution. Bootstrapping can be used to verify how stable the refitting is [@Huang2018]. A more stable refit provides more confidence in the results. It works by making small changes to the mutational profile of a sample. These changes are made by resampling mutations with replacement using the samples own mutational profile as weights. The number of sampled mutations is the same as the number of mutations that was originally in the profile. This process is by default repeated 1000 times. A signature refit is performed for each iteration, resulting in an estimate of the refitting stability. In MutationalPatterns bootstrapping can be done with fit_to_signatures_bootstrapped.

This function can be used with the standard and strict refitting methods described previously. Here we will use the "strict" method on two samples.

contri_boots <- fit_to_signatures_bootstrapped(mut_mat[, c(3, 7)],
  signatures,
  n_boots = 100,
  method = "strict"
)

You can visualize the bootstrapped refitting like this. Each dot is one bootstrap iteration.

plot_bootstrapped_contribution(contri_boots)

You can also visualize this using the relative contribution and a dotplot. Here, the color of the dot shows the percentage of iterations in which the signature is found (contribution > 0), and the size of the dot represents the average contribution of that signature (in the iterations in which the contribution was higher than 0).

plot_bootstrapped_contribution(contri_boots, 
                               mode = "relative", 
                               plot_type = "dotplot")

We can see that SBS1 is relatively stable in the first sample. However, SBS5 is very unstable in the second sample. This instability is likely the result of SBS5 being very flat.

You can also plot the correlation between signatures. A negative correlation between two signatures means that their contributions were high in different bootstrap iterations. Here we will visualize this correlation for one sample.

fig_list <- plot_correlation_bootstrap(contri_boots)
fig_list[[2]]

Here we can see that SBS5 and SBS40 have a negative correlation. This makes sense because they are both flat signatures that are very similar to each other. As a result the refitting process has difficulty distinguishing them.

Similarity between mutational profiles and signatures

Instead of performing NMF or fitting signatures to a profile, you can also look at their similarity. This circumvents the issues that exist with NMF and signature refitting. However, looking at similarities doesn't allow us to separate the different signatures that have contributed to a mutational profile. When multiple signatures have contributed to a profile, the similarities between this profile and the individual signatures can also become diluted.

You can calculate the similarity between two mutational profiles / signatures like this:

cos_sim(mut_mat[, 1], signatures[, 1])

You can also calculate the similarity between multiple mutational profiles / signatures:

cos_sim_samples_signatures <- cos_sim_matrix(mut_mat, signatures)
cos_sim_samples_signatures[1:3, 1:3]

You can visualize this with a heatmap using plot_cosine_heatmap. This function has the same clustering options as plot_contribution_heatmap, which we discussed earlier.

plot_cosine_heatmap(cos_sim_samples_signatures, 
                    cluster_rows = TRUE, cluster_cols = TRUE)

It's also possible to look at the cosine similarities between samples.

cos_sim_samples <- cos_sim_matrix(mut_mat, mut_mat)
plot_cosine_heatmap(cos_sim_samples, cluster_rows = TRUE, cluster_cols = TRUE)

Signature potential damage analysis

Some signatures are more likely than others to have functional effects, by causing "stop gain" or "mismatch" mutations. With MutationalPatterns it's possible to analyze how likely it is for a signature to either cause "stop gain", "mismatch", "synonymous" or "splice site" mutations for a set of genes of interest. Please take into account that this is a relatively basic analysis, that only looks at mutational contexts. Other features like open/closed chromatin are not taken into account. This analysis is meant to give an indication, not a definitive answer, of how damaging a signature might be.

First you need to load a transcription annotation database and make sure some dependencies are installed.

# For example get known genes table from UCSC for hg19 using
# BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
# BiocManager::install("AnnotationDbi")
# BiocManager::install("GenomicFeatures")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

Next, you need to choose a set of genes and create a vector of Entrez gene ids. In this example we used a small set to keep the runtime low, but in practice you can use a larger list of genes, that you are interested in. (The genes used in this example are: P53, KRAS, NRAS, BRAF, BRCA2, CDKN2A, ARID1A, PTEN and TERT.) A useful list of cancer genes can be found here: https://cancer.sanger.ac.uk/cosmic/census.

gene_ids <- c(7157, 3845, 4893, 673, 675, 1029, 8289, 5728, 7015)

Now the ratio of "stop gain", "mismatch", "synonymous" and "splice site" mutations can be determined per genomic context. The total number of possible mutations per context is also given. Finally, a blosum62 score is given for the mismatches. A lower score means that the amino acids in the mismatches are more dissimilar. More dissimilar amino acids are more likely to have a detrimental effect.

contexts <- rownames(mut_mat)
context_mismatches <- context_potential_damage_analysis(contexts, txdb, ref_genome, gene_ids)
head(context_mismatches)

The ratios per context can then be used to get the ratios per signature. Normalized ratios are also given. These were calculated by dividing the ratios in each signature, by the ratios of a completely "flat" signature. A normalized ratio of 2 for "stop gain" mutations, means that a signature is twice as likely to cause "stop gain" mutations, compared to a completely random "flat" signature. The total number of possible mutations per context is multiplied with the signature contribution per context and summed over all contexts. It thus gives a measure of the amount of mutations that a signature could cause.

sig_damage <- signature_potential_damage_analysis(signatures, contexts, context_mismatches)
head(sig_damage)

Using other signature matrixes

So far we have used the SNV signatures from COSMIC. For your convenience we have also included Indel, DBS and transcription strand bias signatures in this package. Additionally, we included signatures from SIGNAL [@Kucab2019, @Degasperi2020]. These signature matrixes can all be loaded using the get_known_signature function. If you use any of these signature matrixes, please cite the associated paper. (The papers are listed in the functions documentation.) A complete list of signature matrixes is shown in the documentation.

You can choose the mutation type like this:

signatures_indel = get_known_signatures(muttype = "indel")
signatures_indel[1:5, 1:5]

It's also possible to include signatures, that might be artifacts. Including these signatures can lead to more overfitting. Therefore we recommend against using them for most analyses. However, these signatures can be usefull to see if your data contains many sequencing artifacts, if you doubt the quality of your data.

signatures_artifacts = get_known_signatures(incl_poss_artifacts = TRUE)
dim(signatures_artifacts)

You can load the SIGNAL reference signatures like this:

signatures_signal = get_known_signatures(source = "SIGNAL")
signatures_signal[1:5, 1:5]

SIGNAL also contains signatures based on drug exposures:

signatures_exposure = get_known_signatures(source = "SIGNAL", sig_type = "exposure")
signatures_exposure[1:5, 1:5]

Finally, SIGNAL contains tissue specific signatures:

signatures_stomach = get_known_signatures(source = "SIGNAL", sig_type = "tissue", tissue_type = "Stomach")
signatures_stomach[1:5, 1:5]

Using an incorrect tissue_type will result in an error. This is usefull, because it shows all possible tissue types. (Not run here, to prevent the error.):

get_known_signatures(source = "SIGNAL", sig_type = "tissue", tissue_type = "?")

The contribution of tissue specific signatures can be converted back to SIGNAL reference signatures. First fit the mutation matrix to tissue specific signatures:

fit_res_tissue <- fit_to_signatures(mut_mat, signatures_stomach)
fit_res_tissue$contribution[1:5, 1:5]

Then convert the contributions to reference signatures:

fit_res_tissue <- convert_sigs_to_ref(fit_res_tissue)
fit_res_tissue$contribution[1:5, 1:5]

Instead of using a signature matrix included in this package, you can also download your own signature matrixes. If you do this you have to make sure that the order of the mutation types is the same as the MutationalPatterns standard. (You can use the match function for this.)

Strand bias analyses

Transcriptional strand bias analysis

For the mutations within genes it can be determined whether the mutation is on the transcribed or non-transcribed strand, which can be used to evaluate the involvement of transcription-coupled repair. To this end, it is determined whether the "C" or "T" base (since by convention we regard base substitutions as C>X or T>X) are on the same strand as the gene definition. Base substitutions on the same strand as the gene definitions are considered "untranscribed", and on the opposite strand of gene bodies as "transcribed", since the gene definitions report the coding or sense strand, which is untranscribed. No strand information is reported for base substitution that overlap with more than one gene body on different strands.

Gene definitions

Start by getting gene definitions for your reference genome:

genes_hg19 <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
genes_hg19

Strand bias profile

You can get transcriptional strand information for all positions in the first VCF object with mut_strand. This function returns "-" for positions outside gene bodies, and positions that overlap with more than one gene on different strands.

strand <- mut_strand(grl[[1]], genes_hg19)
head(strand, 10)

You can make a mutation count matrix with transcriptional strand information (96 trinucleotides * 2 strands = 192 features). NB: only those mutations that are located within gene bodies are counted.

mut_mat_s <- mut_matrix_stranded(grl, ref_genome, genes_hg19)
mut_mat_s[1:5, 1:5]

You can visualize samples from this matrix like this:

plot_192_profile(mut_mat_s[, 1:2])

Strand bias test

You can count the number of mutations on each strand, per tissue, per mutation type:

strand_counts <- strand_occurrences(mut_mat_s, by = tissue)
head(strand_counts)

Next, you can use these counts to perform a Poisson test for strand asymmetry. Multiple testing correction is also performed.

strand_bias <- strand_bias_test(strand_counts)
head(strand_bias)

Plot the mutation spectrum with strand distinction:

ps1 <- plot_strand(strand_counts, mode = "relative")

Plot the effect size (log2(untranscribed/transcribed) of the strand bias. Asteriks indicate significant strand bias. Here we use p-values to plot asterisks. By default fdr is used.

ps2 <- plot_strand_bias(strand_bias, sig_type = "p")

Finally, combine the plots into one figure:

grid.arrange(ps1, ps2)

You can change the significance cutoffs for the fdr and p-values. You can use up to three cutoff levels for each, which changes the number of asteriks in the significant and significant_fdr columns. These asteriks will be used in the plot.

strand_bias_notstrict <- strand_bias_test(strand_counts,
  p_cutoffs = c(0.5, 0.1, 0.05),
  fdr_cutoffs = 0.5
)
plot_strand_bias(strand_bias_notstrict, sig_type = "p")

Replicative strand bias analysis

The involvement of replication-associated mechanisms can be evaluated by testing for a mutational bias between the leading and lagging strand. The replication strand is dependent on the locations of replication origins from which DNA replication is fired. However, replication timing is dynamic and cell-type specific, which makes replication strand determination less straightforward than transcriptional strand bias analysis. Replication timing profiles can be generated with Repli-Seq experiments. Once the replication direction is defined, a strand asymmetry analysis can be performed similarly as the transcription strand bias analysis. The only difference is that you need to use the replication mode for the mut_strand and mut_strand_matrix functions.

Define replication direction

Here we read in an example bed file provided with the package containing replication direction annotation:

repli_file <- system.file("extdata/ReplicationDirectionRegions.bed",
  package = "MutationalPatterns"
)
repli_strand <- read.table(repli_file, header = TRUE)
# Store in GRanges object
repli_strand_granges <- GRanges(
  seqnames = repli_strand$Chr,
  ranges = IRanges(
    start = repli_strand$Start + 1,
    end = repli_strand$Stop
  ),
  strand_info = repli_strand$Class
)
# UCSC seqlevelsstyle
seqlevelsStyle(repli_strand_granges) <- "UCSC"
repli_strand_granges

This GRanges object should have a strand_info metadata column, which contains only two different annotations, e.g. "left" and "right", or "leading" and "lagging". The genomic ranges cannot overlap, to allow only one annotation per location.

The levels of the strand_info metadata in the GRanges object determines the order in which the strands are reported in the mutation matrix that is returned by mut_matrix_stranded, so if you want to count right before left, you can specify this, before you run mut_matrix_stranded:

repli_strand_granges$strand_info <- factor(repli_strand_granges$strand_info,
  levels = c("right", "left")
)

Replication bias analysis

Now that we defined the replication direction, the rest of the analysis is similar to the transcription bias analysis:

You can calculate the strand matrix, counts and bias like this:

mut_mat_s_rep <- mut_matrix_stranded(grl, ref_genome, repli_strand_granges,
  mode = "replication"
)
strand_counts_rep <- strand_occurrences(mut_mat_s_rep, by = tissue)
strand_bias_rep <- strand_bias_test(strand_counts_rep)

And then visualize them:

ps1 <- plot_strand(strand_counts_rep, mode = "relative")
ps2 <- plot_strand_bias(strand_bias_rep)
grid.arrange(ps1, ps2)

Signatures with strand bias

Strand bias can be included in signature analyses. You can for example perform NMF on a mutation count matrix with strand features:

nmf_res_strand <- extract_signatures(mut_mat_s, rank = 2, single_core = TRUE)
colnames(nmf_res_strand$signatures) <- c("Signature A", "Signature B")

Genomic distribution

Mutations are not randomly distributed throughout the genome. With MutationalPatterns you can visualize how mutations are distributed throughout the genome. You can also look at specific genomic regions, such as promoters, CTCF binding sites and transcription factor binding sites. Within these regions you can look for enrichment/depletion of mutations and you can look for differences in the mutational spectra between them.

Rainfall plot

A rainfall plot visualizes mutation types and intermutation distance. Rainfall plots can be used to visualize the distribution of mutations along the genome or a subset of chromosomes. The y-axis corresponds to the distance of a mutation with the previous mutation and is log10 transformed. Drop-downs from the plots indicate clusters or "hotspots" of mutations.

In this example we make a rainfall plot over the autosomal chromosomes for 1 sample:

# Define autosomal chromosomes
chromosomes <- seqnames(get(ref_genome))[1:22]

# Make a rainfall plot
plot_rainfall(grl[[1]],
  title = names(grl[1]),
  chromosomes = chromosomes, cex = 1.5, ylim = 1e+09
)

Define genomic regions

To look at specific types of genomic regions you first need to define them in a named GRangesList. You can use your own genomic region definitions (based on e.g. ChipSeq experiments) or you can use publicly available genomic annotation data, like in the example below.

The following example displays how to download promoter, CTCF binding sites and transcription factor binding sites regions for genome build hg19 from Ensembl using Biocpkg("biomaRt"). For other datasets, see the r Biocpkg("biomaRt") documentation [@Durinck2005]. (Remember to install this package before trying to use it.)

Load the Biocpkg("biomaRt") package.

library(biomaRt)

Download genomic regions. NB: Here we take some shortcuts by loading the results from our example data. The corresponding code for downloading this data can be found above the command we run:

# regulatory <- useEnsembl(biomart="regulation",
#                          dataset="hsapiens_regulatory_feature",
#                          GRCh = 37)

## Download the regulatory CTCF binding sites and convert them to
## a GRanges object.
# CTCF <- getBM(attributes = c('chromosome_name',
#                             'chromosome_start',
#                             'chromosome_end',
#                             'feature_type_name'),
#              filters = "regulatory_feature_type_name",
#              values = "CTCF Binding Site",
#              mart = regulatory)
#
# CTCF_g <- reduce(GRanges(CTCF$chromosome_name,
#                 IRanges(CTCF$chromosome_start,
#                 CTCF$chromosome_end)))

CTCF_g <- readRDS(system.file("states/CTCF_g_data.rds",
  package = "MutationalPatterns"
))

## Download the promoter regions and convert them to a GRanges object.

# promoter = getBM(attributes = c('chromosome_name', 'chromosome_start',
#                                 'chromosome_end', 'feature_type_name'),
#                  filters = "regulatory_feature_type_name",
#                  values = "Promoter",
#                  mart = regulatory)
# promoter_g = reduce(GRanges(promoter$chromosome_name,
#                     IRanges(promoter$chromosome_start,
#                             promoter$chromosome_end)))

promoter_g <- readRDS(system.file("states/promoter_g_data.rds",
  package = "MutationalPatterns"
))

## Download the promoter flanking regions and convert them to a GRanges object.

# flanking = getBM(attributes = c('chromosome_name',
#                                 'chromosome_start',
#                                 'chromosome_end',
#                                 'feature_type_name'),
#                  filters = "regulatory_feature_type_name",
#                  values = "Promoter Flanking Region",
#                  mart = regulatory)
# flanking_g = reduce(GRanges(
#                        flanking$chromosome_name,
#                        IRanges(flanking$chromosome_start,
#                        flanking$chromosome_end)))

flanking_g <- readRDS(system.file("states/promoter_flanking_g_data.rds",
  package = "MutationalPatterns"
))

Combine all genomic regions (GRanges objects) in a named GrangesList:

regions <- GRangesList(promoter_g, flanking_g, CTCF_g)

names(regions) <- c("Promoter", "Promoter flanking", "CTCF")

Make sure that these regions use the same chromosome naming convention as the mutation data:

seqlevelsStyle(regions) <- "UCSC"

Enrichment or depletion of mutations in genomic regions

It is necessary to include a list with GRanges of regions that were surveyed in your analysis for each sample, that is: positions in the genome at which you have enough high quality reads to call a mutation. This can be determined using e.g. CallableLoci by GATK. If you would not include the surveyed area in your analysis, you might for example see a depletion of mutations in a certain genomic region that is solely a result from a low coverage in that region, and therefore does not represent an actual depletion of mutations.

We provided an example surveyed region data file with the package. For simplicity, here we use the same surveyed file for each sample. For a proper analysis, determine the surveyed area per sample and use these in your analysis.

Load the example surveyed region data:

## Get the filename with surveyed/callable regions
surveyed_file <- system.file("extdata/callableloci-sample.bed",
  package = "MutationalPatterns"
)

## Import the file using rtracklayer and use the UCSC naming standard
library(rtracklayer)
surveyed <- import(surveyed_file)
seqlevelsStyle(surveyed) <- "UCSC"

## For this example we use the same surveyed file for each sample.
surveyed_list <- rep(list(surveyed), 9)

First you need to calculate the number of observed and the number of expected mutations in each genomic region for each sample.

distr <- genomic_distribution(grl, surveyed_list, regions)

Next, you can test for enrichment or depletion of mutations in the defined genomic regions using a two-sided binomial test. For this test, the chance of observing a mutation is calculated as the total number of mutations, divided by the total number of surveyed bases. Multiple testing correction is also performed. The significance cutoffs for the fdr and p-values can be changed in the same way as for strand_bias_test. In this example we perform the enrichment/depletion test by tissue type.

distr_test <- enrichment_depletion_test(distr, by = tissue)
head(distr_test)

Finally, you can plot the results. Asteriks indicate significant enrichment/depletion. Here we use p-values to plot asterisks. By default fdr is used.

plot_enrichment_depletion(distr_test, sig_type = "p")

Mutational patterns of genomic regions

Split mutations based on genomic regions

You can also look at the mutational patterns of genomic regions. However, keep in mind that regions with very few mutations will lead to less reliable results.

First you can split the GRangesList containing the mutations based on the defined genomic regions.

grl_region <- split_muts_region(grl, regions)
names(grl_region)

You could now treat these sample/region combinations as completely separate samples. You could for example perform NMF on these, to try to identify signatures that are specific to certain genomic regions.

mut_mat_region <- mut_matrix(grl_region, ref_genome)
nmf_res_region <- extract_signatures(mut_mat_region, rank = 2, nrun = 10, single_core = TRUE)
nmf_res_region <- rename_nmf_signatures(nmf_res_region, 
                                        signatures, 
                                        cutoff = 0.85)
plot_contribution_heatmap(nmf_res_region$contribution, 
                          cluster_samples = TRUE, 
                          cluster_sigs = TRUE)

In this case there don't seem to be any region specific signatures.

Mutation Spectrum

Instead of treating the sample/region combinations as separate samples, you can also plot the spectra per genomic region using the plot_spectrum_region function. The arguments of plot_spectrum can also be used with this function. By default the y-axis shows the number of variants divided by the total number of variants in that sample and genomic region. This way the spectra of regions with very few mutations can be more easily compared to regions with many mutations.

type_occurrences_region <- mut_type_occurrences(grl_region, ref_genome)
plot_spectrum_region(type_occurrences_region)

You can also plot the number of variants divided by the total number of variants in that sample on the y-axis. In this case you don't normalize for the number of variants per genomic region. As you can see below the vast majority of mutations in this example occurred in the "other" region.

plot_spectrum_region(type_occurrences_region, mode = "relative_sample")

Mutational profiles

In addition to plotting the spectra you can also plot a mutational profile. To do this you first need to make a "long" mutation matrix. In this matrix the different genomic regions are considered as different mutational types, instead of as different samples like before.

mut_mat_region <- mut_matrix(grl_region, ref_genome)
mut_mat_long <- lengthen_mut_matrix(mut_mat_region)
mut_mat_long[1:5, 1:5]

You can now plot this using plot_profile_region. The arguments of plot_96_profile can also be used with this function. The options for the y-axis are the same as for plot_spectrum_region. However, by default no normalization is performed for the number of variants per genomic region, because of the often limited number of mutations per mutation type.

plot_profile_region(mut_mat_long[, c(1, 4, 7)])

NB: Since the "mut_mat_long" is a mutation matrix, you could perform NMF on it. This would result in signatures, which will contain different mutation types in different genomic regions.

Mutation density

In the examples above we used known features like promoters for the regions. It's also possible to define regions based on mutation density. You can divide the genome into 3 bins with different mutation density like this:

regions_density <- bin_mutation_density(grl, ref_genome, nrbins = 3)
names(regions_density) <- c("Low", "Medium", "High")

These regions can then be used in the same way as the previous regions. This can be usefull to, for example, compare the spectrum of regions with kataegis with that of the rest of the genome.

grl_region_dens <- split_muts_region(grl, regions_density, include_other = FALSE)

Lesion segregation

Large Watson versus Crick strand asymmetries can sometimes be observed in mutation spectra [@Aitken2020]. This can be the result of many DNA lesions occurring during a single cell cycle. For example, many C>T lesions could occur. If these lesions aren't properly repaired before the next genome duplication, then the resulting sister chromatids will contain the incorrect "T" nucleotides only on their parental strand. Incorrect "A" nucleotides will be incorporated on the newly synthesized strands. These sister chromatids will segregate into different daughter cells, which will have the C>T variants on different strands. The majority of mutations will be either on the Watson or the Crick strand. This process is known as lesion segregation [@Aitken2020].

Healthy human cells are 2n. Therefore, a daughter cell could inherit one copy of a specific chromosome with mutations on the Watson strand and one copy with mutations on the Crick strand. These will cancel each other out and no strand bias will be visible. Because the chromosomes segregate independently from each other, lesion segregation is expected to follow a mendelian inheritance pattern of 1:2:1. 25% of the chromosomes will have mutations on the Watson strand, 25% will have mutations on the Crick strand and 50% will show no Watson versus Crick bias.

Visualizing lesion segregation

You can visualize possible lesion segregation for a single sample. If lesion segregation is present, then it will generally be quite clear that the mutations are not randomly distributed over the strands. In this example no lesion segregation is present. ("+" and "-" are used instead of "Watson" and "Crick" to save space.)

plot_lesion_segregation(grl[[1]])

Calculating lesion segregation

You can also calculate whether lesion segregation is present instead of visualizing it.

You can calculate whether lesion segregation is present using calculate_lesion_segregation. This has the benefit that we can quantify the amount of lesion segregation and generate p-values. However, the generated p-values aren't always 100% reliable as will be discussed below. Therefore, we recommend you to always confirm any suspected lesion segregation by visualizing it.

calculate_lesion_segregation has three different modes. The first mode is based on how often two suqsequent mutations occur on the same strand. The function assumes that when no lesion segregation is present, there is a 50% chance of two subsequent mutations occurring on the same strand. A two-sided binomial test is used to calculate whether the strand between subsequent mutations is switched more often than that. Multiple testing correction is also performed.

lesion_segretation <- calculate_lesion_segregation(grl, sample_names)
head(lesion_segretation)

This statistical test can be influenced by events such as kataegis and local strand assymetries like replication-associated strand bias. As a result the p-value can incorrectly suggest that lesion segregation is present. Therefore, it can be usefull to also look at the fraction of strand switches. In samples with lesion segregation this is generally below 0.4.

By default this mode calculates the lesion segregation for all mutations together. However, a mutational process might cause multiple types of base substitutions, which aren't necessarily considered to be on the same strand. Therefore, it might be usefull to calculate the number of strand switches per mutation type and then sum up the results. In this case the reference genome also needs to be set. We recommend using this when you have a sample with suspected lesion segregation and multiple common types of base substitutions.

lesion_seg_type <- calculate_lesion_segregation(grl,
                                                sample_names,
                                                split_by_type = TRUE,
                                                ref_genome = ref_genome)

The second mode of calculate_lesion_segregation uses the Wald-Wolfowitz test, which was used by @Aitken2020 This test checks whether the Watson and Crick strands are randomly distributed. It's results should generally be similar to the first mode.

lesion_segretation_wald <- calculate_lesion_segregation(grl, sample_names, 
                                                        test = "wald-wolfowitz")
head(lesion_segretation_wald)

This statistical test can also be influenced by events such as kataegis and local strand assymetries like replication-associated strand bias.

The third mode of calculate_lesion_segregation can calculate the rl20 value and the associated genomic span, which together are somewhat less sensitive to events like kataegis.

A rl20 value of 6 means that at least 20% of mutations are in a strand specific run of 6 or more consecutive mutations. The genomic span is the part of the genome covered by these runs. If the rl20 is high and a decent part of the genome is covered by the strand specific runs, then this provides strong evidence of lesion segregation. A high rl20, combined with a low genomic span (<5%) is indicative of local clustering events like kataegis [@Aitken2020]. A downside of this method is that it doesn't generate a p-value. In general, we recommend you to use the first or second mode of calculate_lesion_segregation to get a p-value and the third mode to check if you are looking at a genome wide process or a local process like kataegis.

lesion_segretation_rl20 <- calculate_lesion_segregation(grl,
  sample_names,
  test = "rl20",
  ref_genome = ref_genome,
  chromosomes = chromosomes
)
head(lesion_segretation_rl20)

A note on the graphics

The plots made with this package are all made using r CRANpkg("ggplot2") [@Wickham2016]. This means that all the plots (except for the plots with dendograms) are highly customizable. You can for example change the size and text orientation of the y-axis.

p <- plot_spectrum(type_occurrences, legend = FALSE)
p_axis <- p +
  theme(axis.text.y = element_text(size = 14, angle = 90))

You can also change the entire theme of the plot.

p_theme <- p +
  theme_classic()

grid.arrange(p, p_axis, p_theme, ncol = 3, widths = c(3, 3, 3))

More information on r CRANpkg("ggplot2") is available here. A list of themes is available here.

Session Info

sessionInfo()

References

MutationalPatterns documentation built on Nov. 14, 2020, 2:03 a.m.