---
title: "Mutational Signature Comprehensive Analysis Toolkit"
author:
- name: Aaron Chevalier
  affiliation: &id Boston University School of Medicine
- name: Joshua Campbell
  affiliation: *id
  email: camp@bu.edu
output:
  BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{Mutational Signature Comprehensive Analysis Toolkit}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE, results = "asis"}
BiocStyle::markdown()
knitr::opts_chunk$set(echo = TRUE, dev = "png")
```

# Vignettes
Due to the range of features that are available in musicatk, we have prepared several vignettes to help you get started, which are all available at [https://www.camplab.net/musicatk/](https://www.camplab.net/musicatk/)

#Introduction
A variety of exogenous exposures or endogenous biological processes can contribute to the overall mutational load observed in human tumors. Many different mutational patterns, or “mutational signatures”, have been identified across different tumor types. These signatures can provide a record of environmental exposure and can give clues about the etiology of carcinogenesis. The Mutational Signature Comprehensive Analysis Toolkit (musicatk) has utilities for extracting variants from a variety of file formats, contains multiple methods for discovery of novel signatures or prediction of known signatures, as well as many types of downstream visualizations for exploratory analysis. This package has the ability to parse and combine multiple motif classes in the mutational signature discovery or prediction processes. Mutation motifs include single base substitutions (SBS), double base substitutions (DBS), insertions (INS) and deletions (DEL).

# Installation
Currently musicatk can be installed from on Bioconductor using the following code:

```{r, eval= FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("musicatk")
```

To install the latest version from Github, use the following code:

```{r, eval = FALSE}
if (!requireNamespace("devtools", quietly = TRUE)) {
  install.packages("devtools")
}

library(devtools)
install_github("campbio/musicatk")
```

The package can be loaded using the `library` command.

```{r, eval = TRUE, message = FALSE}
library(musicatk)
```

# Setting up a musica object
In order to discover or predict mutational signatures, we must first set up
our musica object by 1) extracting variants from files or objects such as
VCFs and MAFs, 2) selecting the appropriate reference genome 3) creating a 
musica object, and 4) building a count tables for our variants of interest.
Alternatively, a musica object can be created directly from a count table.

## Extracting variants
Variants can be extracted from various formats using the following functions:

* The `extract_variants_from_vcf_file()` function will extract variants from a [VCF](https://samtools.github.io/hts-specs/) file. The file will be imported using the readVcf function from the [VariantAnnotation](https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html) package and then the variant information will be extracted from this object.
* The `extract_variants_from_vcf()` function extracts variants from a `CollapsedVCF` or `ExpandedVCF` object from the [VariantAnnotation](https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html) package.
* The `extract_variants_from_maf_file()` function will extract variants from a file in [Mutation Annotation Format (MAF)](https://docs.gdc.cancer.gov/Data/File_Formats/MAF_Format/) used by TCGA.
* The `extract_variants_from_maf()` function will extract variants from a MAF object created by the [maftools](https://www.bioconductor.org/packages/release/bioc/html/maftools.html) package. 
* The `extract_variants_from_matrix()` function will get the information from a matrix or data.frame like object that has columns for the chromosome, start position, end position, reference allele,  mutation allele, and sample name. 
* The `extract_variants()` function will extract variants from a list of objects. These objects can be any combination of VCF files, VariantAnnotation objects, MAF files, MAF objects, and data.frame objects.

Below are some examples of extracting variants from MAF and VCF files:

```{r extract_variants, message = FALSE}
# Extract variants from a MAF File
lusc_maf <- system.file("extdata", "public_TCGA.LUSC.maf", package = "musicatk")
lusc.variants <- extract_variants_from_maf_file(maf_file = lusc_maf)

# Extract variants from an individual VCF file
luad_vcf <- system.file("extdata", "public_LUAD_TCGA-97-7938.vcf",
  package = "musicatk"
)
luad.variants <- extract_variants_from_vcf_file(vcf_file = luad_vcf)

# Extract variants from multiple files and/or objects
melanoma_vcfs <- list.files(system.file("extdata", package = "musicatk"),
  pattern = glob2rx("*SKCM*vcf"), full.names = TRUE
)
variants <- extract_variants(c(lusc_maf, luad_vcf, melanoma_vcfs))
```

## Choosing a genome
musicatk uses [BSgenome](https://bioconductor.org/packages/release/bioc/html/BSgenome.html) objects to access genome sequence information that flanks each mutation which is used bases for generating mutation count tables. BSgenome objects store full genome sequences for different organisms. A full list of supported organisms can be obtained by running `available.genomes()` after loading the BSgenome library. Custom genomes can be forged as well (see [BSgenome](https://bioconductor.org/packages/release/bioc/html/BSgenome.html) documentation). musicatk provides a utility function called `select_genome()` to allow users to quickly select human genome build versions "hg19" and "hg38" or mouse genome builds "mm9" and "mm10". The reference sequences for these genomes are in UCSC format (e.g. chr1).

```{r select_genome}
g <- select_genome("hg38")
```

## Creating a musica object
The last preprocessing step is to create an object with the variants and the genome using the `create_musica_from_variants` function. This function will perform checks to ensure that the chromosome names and reference alleles in the input variant object match those in supplied BSgenome object. These checks can be turned off by setting `check_ref_chromosomes = FALSE` and `check_ref_bases = FALSE`, respectively. This function also looks for adjacent single base substitutions (SBSs) and will convert them to double base substitutions (DBSs). To disable this automatic conversion, set `convert_dbs = FALSE`. 

```{r create_musica}
musica <- create_musica_from_variants(x = variants, genome = g)
```


# Creating mutation count tables
Motifs are the building blocks of mutational signatures. Motifs themselves are 
a mutation combined with other genomic information. For instance, **SBS96** 
motifs are constructed from an SBS mutation and one upstream and one downstream 
base sandwiched together. We build tables by counting these motifs for each 
sample. 
```{r build_tables}
build_standard_table(musica, g = g, table_name = "SBS96")
```

Here is a list of mutation tables that can be created by setting the 
`table_name` parameter in the `build_standard_table` function:

* SBS96 - Motifs are the six possible single base pair mutation types times the 
four possibilities each for upstream and downstream context bases (4*6*4 = 96 
motifs)
* SBS192_Trans - Motifs are an extension of SBS96 multiplied by the 
transcriptional strand (translated/untranslated), can be specified with 
`"Transcript_Strand"`.
* SBS192_Rep - Motifs are an extension of SBS96 multiplied by the 
replication strand (leading/lagging), can be specified with 
`"Replication_Strand"`.
* DBS - Motifs are the 78 possible double-base-pair substitutions 
* INDEL - Motifs are 83 categories intended to capture different categories of 
indels based on base-pair change, repeats, or microhomology, insertion or 
deletion, and length.

```{r combine_tables}
data(dbs_musica)
build_standard_table(dbs_musica, g, "SBS96", overwrite = TRUE)
build_standard_table(dbs_musica, g, "DBS78", overwrite = TRUE)

# Subset SBS table to DBS samples so they cam be combined
count_tables <- tables(dbs_musica)
overlap_samples <- which(colnames(count_tables$SBS96@count_table) %in%
                           colnames(count_tables$DBS78@count_table))
count_tables$SBS96@count_table <-
  count_tables$SBS96@count_table[, overlap_samples]
tables(dbs_musica) <- count_tables

combine_count_tables(
  musica = dbs_musica, to_comb = c("SBS96", "DBS78"),
  name = "sbs_dbs", description = "An example combined
                     table, combining SBS96 and DBS", overwrite = TRUE
)
```


```{r}
annotate_transcript_strand(musica, "19", build_table = FALSE)
build_custom_table(
  musica = musica, variant_annotation = "Transcript_Strand",
  name = "Transcript_Strand",
  description = "A table of transcript strand of variants",
  data_factor = c("T", "U"), overwrite = TRUE
)
```


Different count tables can be combined into one using the `combine_count_tables`
function. For example, the SBS96 and the DBS tables could be combined and 
mutational signature discovery could be performed across both mutations 
modalities. Tables with information about the same types of variants (e.g. 
two related SBS tables) should generally not be combined and used together.

Custom count tables can be created from user-defined mutation-level annotations
using the `build_custom_table` function.

# Creating a musica object directly from a count table
If a count table is already available, a musica object can be created directly
without need for a variant file and building tables.

```{r create_musica_from_counts}
luad_count_table_path <- system.file("extdata", "luad_tcga_count_table.csv",
  package = "musicatk"
)
luad_count_table <- as.matrix(read.csv(luad_count_table_path))

musica_from_counts <- create_musica_from_counts(luad_count_table, "SBS96")
```

# Discovering Signatures and Exposures
Mutational signature discovery is the process of deconvoluting a matrix
containing counts for each mutation type in each sample into two matrices: 1)
a **signature** matrix containing the probability of each mutation motif in
signature and 2) an **exposure** matrix containing the estimated counts of each
signature in each sample. Discovery and prediction results are saved in the 
result_list slot of a `musica` object.

```{r discover_sigs}
discover_signatures(
  musica = musica, modality = "SBS96", num_signatures = 3,
  algorithm = "nmf", model_id = "ex_result", nstart = 10
)
```

Supported signature discovery algorithms include:

* Non-negative matrix factorization (nmf)
* Latent Dirichlet Allocation (lda)

Both have built-in `seed` capabilities for reproducible results, `nstarts` for 
multiple independent chains from which the best final result will be chosen. 
NMF also allows for parallel processing via `par_cores`.

To get the signatures or exposures from a result within a musica object, the 
following functions can be used:

```{r result_accessors}
# Extract the exposure matrix
expos <- exposures(musica, "result", "SBS96", "ex_result")
expos[1:3, 1:3]

# Extract the signature matrix
sigs <- signatures(musica, "result", "SBS96", "ex_result")
sigs[1:3, 1:3]
```

## Determining an appropriate k value
The k value is the number of signatures that are predicted from a discovery
method. To help determine an appropriate k value, the `compare_k_vals` function
can be used to compare to the stability and error associated with various k
values. Generally, 100 replicates is suggested, as well as a larger span of k
values to test. Here, fewer k values and replicates are used for simplicity.

```{r compare_k_vals}
k_comparison <- compare_k_vals(musica, "SBS96",
  reps = 5, min_k = 2, max_k = 4,
  algorithm = "nmf"
)
```

# Plotting

## Signatures
Barplots showing the probability of each mutation type in each signature can 
be plotted with the `plot_signatures` function:

```{r, plot_sigs}
plot_signatures(musica, model_id = "ex_result", modality = "SBS96")
```

By default, the scales on the y-axis are forced to be the same across all 
signatures. This behavior can be turned off by setting `same_scale = FALSE`. 
Signatures can be re-named based on prior knowledge and displayed in the plots:

```{r, name_sigs}
name_signatures(musica, "ex_result", c("Smoking", "APOBEC", "UV"))
plot_signatures(musica, model_id = "ex_result", modality = "SBS96")
```


## Exposures

Barplots showing the exposures in each sample can be plotted with the 
`plot_exposures` function:

```{r exposures_raw}
plot_exposures(musica, "ex_result", plot_type = "bar")
```

The proportion of each exposure in each tumor can be shown by setting `proportional = TRUE`:
```{r exposures_prop}
plot_exposures(musica, "ex_result", plot_type = "bar", proportional = TRUE)
```

Counts for individual samples can also be plotted with the `plot_sample_counts` function:

```{r sample_counts}
samples <- sample_names(musica)
plot_sample_counts(musica, sample_names = samples[c(3, 4, 5)],
                   modality = "SBS96")
```

## Comparison to external signatures (e.g. COSMIC)
A common analysis is to compare the signatures estimated in a dataset to those generated in other datasets or to those in the [COSMIC database](https://cancer.sanger.ac.uk/cosmic/signatures). We have a set of functions that can be used to easily perform pairwise correlations between signatures. The `compare_results` functions compares the signatures between two models in the same or different `musica` objects. The `compare_cosmic_v2` will correlate the signatures between a model and the SBS signatures in COSMIC V2. For example:

```{r compare_cosmic}
compare_cosmic_v2(musica, "ex_result", threshold = 0.75)
```

In this example, our Signatures 1 and 2 were most highly correlated to COSMIC Signature 7 and 4, respectively, so this may indicate that samples in our dataset were exposed to UV radiation or cigarette smoke. Only pairs of signatures who have a correlation above the `threshold` parameter will be returned. If no pairs of signatures are found, then you may want to consider lowering the threshold. The default correlation metric is the cosine similarity, but this can be changed to using 1 - Jensen-Shannon Divergence by setting 
`metric = "jsd"`  Signatures can also be correlated to those in the COSMIC V3 database using the `compare_cosmic_v3` function.

# Predicting exposures using pre-existing signatures

Instead of discovering mutational signatures and exposures from a dataset *de novo*, one might get better results by predicting the exposures of signatures that have been previously estimated in other datasets. We incorporate several methods for estimating exposures given a set of pre-existing signatures. For example, we can predict exposures for COSMIC signatures 1, 4, 7, and 13 in our current dataset:

```{r predict_cosmic}
# Load COSMIC V2 data
data("cosmic_v2_sigs")

# Predict pre-existing exposures using the "lda" method
predict_exposure(
  musica = musica, modality = "SBS96",
  signature_res = cosmic_v2_sigs, model_id = "ex_pred_result",
  signatures_to_use = c(1, 4, 7, 13), algorithm = "lda"
)

# Plot exposures
plot_exposures(musica, "ex_pred_result", plot_type = "bar")
```

The `cosmic_v2_sigs` object is a `result_model` object containing COSMIC V2 signatures without any sample or exposure information. Note that if `signatures_to_use` is not supplied by the user, then exposures for all signatures in the result object will be estimated. We can predict exposures for samples in any `musica` object from any `result_model` object as long as the same mutation schema was utilized. 

We can list which signatures are present in each tumor type according to the [COSMIC V2 database](https://cancer.sanger.ac.uk/cosmic/signatures_v2.tt). For example, we can find which signatures are present in lung cancers:

```{r subtype_map}
cosmic_v2_subtype_map("lung")
```

We can predict exposures for samples in a `musica` object using the signatures from any `result_model` object. Just for illustration, we will predict exposures using the estimated signatures from model we created earlier:

```{r predict_previous}
model <- get_model(musica, "result", "SBS96", "ex_result")
predict_exposure(
  musica = musica, modality = "SBS96",
  signature_res = model, algorithm = "lda",
  model_id = "pred_our_sigs"
)
```

Of course, this example is not very useful because we are predicting exposures using signatures that were learned on the same set of samples. Most of the time, you want to give the `signature_res` parameter a `result_model` object that was generated using independent samples from those in the input `musica` object. As mentioned above, different signatures in different models can be compared to each other using the `compare_results` function:

```{r predict_compare}
compare_results(musica,
  model_id = "ex_pred_result",
  other_model_id = "pred_our_sigs", threshold = 0.60
)
```

# Comparing samples between groups using Sample Annotations

## Adding sample annotations

We first must add an annotation to the `musica` object
```{r annotations}
annot <- read.table(system.file("extdata", "sample_annotations.txt",
                                package = "musicatk"),
                    sep = "\t", header = TRUE)
samp_annot(musica, "Tumor_Subtypes") <- annot$Tumor_Subtypes
```

Note that the annotations can also be added directly the `musica` object in the beginning using the same function: `samp_annot(musica, "Tumor_Subtypes") <- annot$Tumor_Subtypes`. These annotations will then automatically be included in any down-stream result object.

* **Be sure that the annotation vector being supplied is in the same order as the samples in the `musica` object.** You can view the sample order with the `sample_names` function:

```{r sample_names}
sample_names(musica)
```

## Plotting exposures by a Sample Annotation

As mentioned previously, the `plot_exposures` function can plot sample exposures in a `musica_result` object. It can group exposures by either a sample annotation or by a signature by setting the `group_by` parameter. Here will will group the exposures by the `Tumor_Subtype` annotation:

```{r plot_exposures_by_subtype}
plot_exposures(musica, "ex_result",
  plot_type = "bar", group_by = "annotation",
  annotation = "Tumor_Subtypes"
)
```

The distribution of exposures with respect to annotation can be viewed using boxplots by setting `plot_type = "box"` and `group_by = "annotation"`:

```{r plot_exposures_box_annot}
plot_exposures(musica, "ex_result", plot_type = "box", group_by = "annotation",
               annotation = "Tumor_Subtypes")
```

Note that the name of the annotation must be supplied via the `annotation` parameter. Boxplots can be converted to violin plots by setting `plot_type = "violin"`. To compare the level of each exposure across sample groups within a signature, we can set `group_by = "signature"` and `color_by = "annotation"`:

```{r plot_exposures_box_sig}
plot_exposures(musica, "ex_result",
  plot_type = "box", group_by = "signature",
  color_by = "annotation", annotation = "Tumor_Subtypes"
)
```

## Visualizing samples in 2D using UMAP

The `create_umap` function embeds samples in 2 dimensions using the `umap` function from the [uwot](https://cran.r-project.org/web/packages/uwot/index.html) package. The major parameters for fine tuning the UMAP are `n_neighbors`, `min_dist`, and `spread`. See `?uwot::umap` for more information on these parameters.

```{r umap_create}
create_umap(musica, "ex_result")
```

The `plot_umap` function will generate a scatter plot of the UMAP coordinates. The points of plot will be colored by the level of a signature by default:

```{r umap_plot}
plot_umap(musica, "ex_result")
```

By default, the exposures for each sample will share the same color scale. However, exposures for some signatures may have really high levels compared to others. To make a plot where exposures for each signature will have their own color scale, you can set `same_scale = FALSE`:

```{r umap_plot_same_scale}
plot_umap(musica, "ex_result", same_scale = FALSE)
```

Lastly, points can be colored by a Sample Annotation by setting `color_by = "annotation"` and `annotation` to the name of the annotation:

```{r umap_plot_annot}
plot_umap(musica, "ex_result",
  color_by = "annotation",
  annotation = "Tumor_Subtypes", add_annotation_labels = TRUE
)
```

When `add_annotation_labels = TRUE`, the centroid of each group is identified 
using medians and the labels are plotted on top of the centroid.

# Use of Plotly in plotting

plot_signatures, plot_exposures, and plot_umap, all have builty in ggplotly 
capabilities. Simply specifying `plotly = TRUE` enables interactive plots 
that allows examination of individuals sections, zooming and resizing, and 
turning on and off annotation types and legend values.

```{r plotly}
plot_signatures(musica, "ex_result", plotly = TRUE)
plot_exposures(musica, "ex_result", plotly = TRUE)
plot_umap(musica, "ex_result", plotly = TRUE)
```

# Note on reproducibility

Several functions make use of stochastic algorithms or procedures which 
require the use of random number generator (RNG) for simulation or sampling. 
To maintain reproducibility, all these functions should be called using 
```set_seed(1)``` or ```withr::with_seed(seed, function())``` to make sure 
same results are generatedeach time one of these functions is called. Using
with_seed for reproducibility has the advantage of not interfering with any 
other user seeds, but using one or the other is important for several functions 
including *discover_signatures*, *predict_exposure*, and *create_umap*, as 
these functions use stochastic processes that may produce different results 
when run multiple times with the same settings.

```{r reproducible_prediction}
seed <- 1
withr::with_seed(seed, predict_exposure(
  musica = musica, modality = "SBS96",
  signature_res = model, algorithm = "lda",
  model_id = "reproducible_ex"
))
```


# Session Information
```{r session}
sessionInfo()
```