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)

##  [1] "A>C" "A>G" "C>T" "A>G" "G>T" "T>A" "T>C" "G>A" "G>A" "C>A" "G>A" "G>T"

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)

##  [1] "T>G" "T>C" "C>T" "T>C" "C>A" "T>A" "T>C" "C>T" "C>T" "C>A" "C>T" "C>A"

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)

##  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr2  chr2  chr2 
## "CAG" "AAC" "ACA" "AAG" "TGA" "GTT" "ATT" "CGC" "AGC" "ACA" "CGT" "GGA"

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)

## $types
##  [1] "T>G" "T>C" "C>T" "T>C" "C>A" "T>A" "T>C" "C>T" "C>T" "C>A" "C>T" "C>A"
## 
## $context
##  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr1  chr2  chr2  chr2 
## "CTG" "GTT" "ACA" "CTT" "TCA" "GTT" "ATT" "GCG" "GCT" "ACA" "ACG" "TCC"

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

##            C>A C>G C>T T>A T>C T>G C>T at CpG C>T other
## colon1      28   5 109  12  30  12         59        50
## colon2      77  29 345  36  90  21        209       136
## colon3      79  19 243  25  61  23        165        78
## intestine1  19   8  74  19  26   4         33        41
## intestine2 118  49 423  57 126  27        258       165
## intestine3  54  27 298  32  67  22        192       106
## liver1      43  22  94  30  77  34         18        76
## liver2     144  93 274 103 209  73         48       226
## liver3      39  28  61  15  32  23          7        54

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 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 trinucleotide 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)

##         colon1 colon2 colon3 intestine1 intestine2 intestine3 liver1 liver2 liver3
## A[C>A]A      3     10     10          5         19          6      8     10      3
## A[C>A]C      0      3      3          1          8          4      1      8      2
## A[C>A]G      2      3      3          1          4          0      1      6      2
## A[C>A]T      0      2      9          0          9          2      2     12      2
## C[C>A]A      1      8      5          0          8          7      2     15      3
## C[C>A]C      2      5      3          1          3          2      1     15      2

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 useful if you have a large number of mutations.

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

##           colon1 colon2 colon3 intestine1 intestine2 intestine3 liver1 liver2 liver3
## AA[C>A]AA      0      4      1          2          9          4      0      4      0
## AA[C>A]AC      0      0      2          0          1          0      1      0      0
## AA[C>A]AG      0      0      1          0          2          0      1      0      0
## AA[C>A]AT      0      0      1          1          2          0      0      0      0
## AA[C>A]CA      0      0      0          0          0          0      0      2      1
## AA[C>A]CC      0      1      0          0          2          2      0      1      0

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)])

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 NMF package (Gaujoux and Seoighe 2010). 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 pseudocount 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")

## Compute NMF rank= 2  ... + measures ... OK
## Compute NMF rank= 3  ... + measures ... OK
## Compute NMF rank= 4  ... + measures ... OK
## Compute NMF rank= 5  ... + measures ... OK

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 algorithm: 'brunet'

## Multiple runs: 10

## Mode: sequential [sapply]

## Runs: 1* 2 3 4 5* 6 7 8* 9 10 ... DONE
## System time:
##    user  system elapsed 
##   7.318   0.000   7.319

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)

## NMF algorithm: 'brunet'

## Multiple runs: 10

## Mode: sequential [sapply]

## Runs: 1* 2 3 4 5 6 7* 8 9 10 ... DONE
## System time:
##    user  system elapsed 
##   7.289   0.000   7.290

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 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 useful 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.2) (Alexandrov et al. 2020). (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)

## [1] "SBS5-like" "SBS1-like"

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

## [1] "SBS5-like" "SBS1-like"

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

## [1] "intestine1" "intestine3" "colon3"     "intestine2" "colon1"     "colon2"     "liver3"    
## [8] "liver1"     "liver2"

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)

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. It can work with a SNV, indel, DBS or other type of count matrix.

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 Degasperi et al. (2020). 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 Alexandrov et al. (2020). 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 and depends on the data. Here we use a cutoff of 0.004. Decreasing this number will make the refitting less strict, while increasing it will make the refitting more strict. Trying out different values can sometimes be useful to achieve the best results.

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"
)

By default fit_to_signatures_strict uses the “backwards” selection approach described above. However, it is also possible to use a “best subset” approach. The benefit of this method is that it can be more accurate than the “backwards” approach. However, it becomes computationally infeasible when using many signatures. Therefore it should only be used on small signature sets (max 10-15 signatures), like tissue specific signatures. The “best subset” approach works similarly to the “backwards” approach. This approach again starts with a standard refit. The refitting is then repeated for each combination of n-1 signatures, where n is the total number of signatures. In other words, if you started with 10 signatures, the refitting is repeated 10 times, with a different signature being removed each time. The combination of signatures that has the best cosine similarity between the original and reconstructed profile is chosen. This is done in an iterative fashion for n-2, n-3, ect. You stop removing signatures when the difference between two iterations becomes bigger than a certain cutoff, just like with the backwards method.

We randomly selected a few signatures for this example, to keep the runtime low. In practice, signatures should be selected based on prior knowledge.

best_subset_refit <- fit_to_signatures_strict(mut_mat, signatures[,1:5], max_delta = 0.002, method = "best_subset")

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)

## Combined the following two signatures: SBS26, SBS12

## Combined the following two signatures: SBS36, SBS18

## Combined the following two signatures: SBS92, SBS5

## Combined the following two signatures: SBS40, SBS3

## Combined the following two signatures: SBS10d, SBS10a

## Combined the following two signatures: SBS10d;SBS10a, SBS10c

## Combined the following two signatures: SBS15, SBS6

## Combined the following two signatures: SBS29, SBS24

## Combined the following two signatures: SBS94, SBS4

## Combined the following two signatures: SBS23, SBS19

## Combined the following two signatures: SBS26;SBS12, SBS37

## Combined the following two signatures: SBS86, SBS39

## Combined the following two signatures: SBS40;SBS3, SBS92;SBS5

## Combined the following two signatures: SBS94;SBS4, SBS8

## Combined the following two signatures: SBS23;SBS19, SBS31

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 (Huang, Wojtowicz, and Przytycka 2018). 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. (We only use 50 bootstraps here to reduce the run time and figure size.)

contri_boots <- fit_to_signatures_bootstrapped(mut_mat[, c(3, 7)],
  signatures,
  n_boots = 50,
  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.

Gene definitions

Start by getting gene definitions for your reference genome:

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

##   403 genes were dropped because they have exons located on both strands of the same
##   reference sequence or on more than one reference sequence, so cannot be represented
##   by a single genomic range.
##   Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList object, or
##   use suppressMessages() to suppress this message.

genes_hg19

## GRanges object with 23056 ranges and 1 metadata column:
##         seqnames              ranges strand |     gene_id
##            <Rle>           <IRanges>  <Rle> | <character>
##       1    chr19   58858172-58874214      - |           1
##      10     chr8   18248755-18258723      + |          10
##     100    chr20   43248163-43280376      - |         100
##    1000    chr18   25530930-25757445      - |        1000
##   10000     chr1 243651535-244006886      - |       10000
##     ...      ...                 ...    ... .         ...
##    9991     chr9 114979995-115095944      - |        9991
##    9992    chr21   35736323-35743440      + |        9992
##    9993    chr22   19023795-19109967      - |        9993
##    9994     chr6   90539619-90584155      + |        9994
##    9997    chr22   50961997-50964905      - |        9997
##   -------
##   seqinfo: 93 sequences (1 circular) from hg19 genome

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)

##  [1] transcribed   -             transcribed   -             untranscribed -            
##  [7] -             -             untranscribed -            
## Levels: untranscribed transcribed -

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]

##                       colon1 colon2 colon3 intestine1 intestine2
## A[C>A]A-untranscribed      0      0      0          0          4
## A[C>A]A-transcribed        1      1      2          4          3
## A[C>A]C-untranscribed      0      0      1          1          1
## A[C>A]C-transcribed        0      0      0          0          1
## A[C>A]G-untranscribed      1      0      0          0          0

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)

## # A tibble: 6 × 5
##   group type  strand        no_mutations relative_contribution
##   <fct> <chr> <chr>                <dbl>                 <dbl>
## 1 colon C>A   transcribed             32                0.0731
## 2 colon C>A   untranscribed           23                0.0525
## 3 colon C>G   transcribed             11                0.0251
## 4 colon C>G   untranscribed           10                0.0228
## 5 colon C>T   transcribed            134                0.306 
## 6 colon C>T   untranscribed          115                0.263

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)

## # A tibble: 6 × 10
##   group type  transcribed untranscribed total ratio p_poisson significant   fdr significant_fdr
##   <fct> <chr>       <dbl>         <dbl> <dbl> <dbl>     <dbl> <chr>       <dbl> <chr>          
## 1 colon C>A            32            23    55  1.39     0.281 ""          0.692 ""             
## 2 colon C>G            11            10    21  1.1      1     ""          1     ""             
## 3 colon C>T           134           115   249  1.17     0.254 ""          0.692 ""             
## 4 colon T>A            12             9    21  1.33     0.664 ""          0.910 ""             
## 5 colon T>C            36            32    68  1.12     0.716 ""          0.910 ""             
## 6 colon T>G            15             9    24  1.67     0.307 ""          0.692 ""

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")

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

## GRanges object with 1993 ranges and 1 metadata column:
##          seqnames            ranges strand | strand_info
##             <Rle>         <IRanges>  <Rle> | <character>
##      [1]     chr1   2133001-3089000      * |       right
##      [2]     chr1   3089001-3497000      * |        left
##      [3]     chr1   3497001-4722000      * |       right
##      [4]     chr1   5223001-6428000      * |        left
##      [5]     chr1   6428001-7324000      * |       right
##      ...      ...               ...    ... .         ...
##   [1989]     chrY 23997001-24424000      * |       right
##   [1990]     chrY 24424001-28636000      * |        left
##   [1991]     chrY 28636001-28686000      * |       right
##   [1992]     chrY 28686001-28760000      * |        left
##   [1993]     chrY 28760001-28842000      * |       right
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

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)

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)

##  [1] "colon1.Promoter"              "colon1.Promoter flanking"    
##  [3] "colon1.CTCF"                  "colon1.Other"                
##  [5] "colon2.Promoter"              "colon2.Promoter flanking"    
##  [7] "colon2.CTCF"                  "colon2.Other"                
##  [9] "colon3.Promoter"              "colon3.Promoter flanking"    
## [11] "colon3.CTCF"                  "colon3.Other"                
## [13] "intestine1.Promoter"          "intestine1.Promoter flanking"
## [15] "intestine1.CTCF"              "intestine1.Other"            
## [17] "intestine2.Promoter"          "intestine2.Promoter flanking"
## [19] "intestine2.CTCF"              "intestine2.Other"            
## [21] "intestine3.Promoter"          "intestine3.Promoter flanking"
## [23] "intestine3.CTCF"              "intestine3.Other"            
## [25] "liver1.Promoter"              "liver1.Promoter flanking"    
## [27] "liver1.CTCF"                  "liver1.Other"                
## [29] "liver2.Promoter"              "liver2.Promoter flanking"    
## [31] "liver2.CTCF"                  "liver2.Other"                
## [33] "liver3.Promoter"              "liver3.Promoter flanking"    
## [35] "liver3.CTCF"                  "liver3.Other"

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 algorithm: 'brunet'

## Multiple runs: 10

## Mode: sequential [sapply]

## Runs: 1* 2 3* 4 5 6 7 8 9 10 ... DONE
## System time:
##    user  system elapsed 
##   9.355   0.004   9.360

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]

##                  colon1 colon2 colon3 intestine1 intestine2
## A[C>A]A_Promoter      0      0      0          0          1
## A[C>A]C_Promoter      0      0      0          0          0
## A[C>A]G_Promoter      0      0      0          0          0
## A[C>A]T_Promoter      0      0      1          0          0
## C[C>A]A_Promoter      0      0      0          0          0

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 useful 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)