Paths(signatures$SignExposures)
Motivation: Cancer is an evolutionary process driven by continuous acquisition of genetic variations in individual cells. The diversity and complexity of somatic mutational processes is a conspicuous feature orchestrated by DNA damage agents and repair processes, including exogenous or endogenous mutagen exposures, defects in DNA mismatch repair and enzymatic modification of DNA. The identification of the underlying mutational processes is central to understanding of cancer origin and evolution.
The signeR package focuses on the estimation and further analysis of mutational signatures. The functionalities of this package can be divided into three categories. First, it provides tools to process VCF files and generate matrices of SNV mutation counts and mutational opportunities, both defined according to a 3bp context (mutation site and its neighbouring 3' and 5' bases). Second, these count matrices are considered as input for the estimation of the underlying mutational signatures and the number of active mutational processes. Third, the package provides tools to correlate the activities of those signatures with other relevant information such as clinical data, in order to draw conclusions about the analysed genome samples, which can be useful for clinical applications. These include the Differential Exposure Score and the a posteriori sample classification.
Although signeR is intended for the estimation of mutational signatures, it actually provides a full Bayesian treatment to the non-negative matrix factorisation (NMF) model. Further details about the method can be found in Rosales & Drummond et al., 2016 (see section 6.1 below).
This vignette briefly explains the use of signeR through examples.
Before installing, please make sure you have the latest version of R and Bioconductor installed.
To install signeR, start R and enter:
install.packages("BiocManager") BiocManager::install("signeR")
For more information, see this page.
Once installed the library can be loaded as
library(signeR)
signeR takes as input a count matrix of samples x features. Each feature is usually a SNV mutation within a 3bp context (96 features, 6 types of SNV mutations and 4 possibilities for the bases at each side of the SNV change). Optionally, an opportunity matrix can also be provided containing the count frequency of the features in the whole analyzed region for each sample. Although not required, this argument is highly recommended because it allows signeR to normalize the features frequency over the analyzed region.
Input matrices can be read both from a VCF or a tab-delimited files, as described next.
The VCF file format is the most common format for storing genetic variations, the signeR package includes a utility function for generating a count matrix from the VCF:
library(VariantAnnotation) # BSgenome, equivalent to the one used on the variant call library(BSgenome.Hsapiens.UCSC.hg19) vcfobj <- readVcf("/path/to/a/file.vcf", "hg19") mut <- genCountMatrixFromVcf(BSgenome.Hsapiens.UCSC.hg19, vcfobj)
This function will generate a matrix of mutation counts for each sample in the provided VCF.
If you have one vcf per sample you can combine the results into a single matrix like this:
mut = matrix(ncol=96,nrow=0) for(i in vcf_files) { vo = readVcf(i, "hg19") # sample name (should pick up from the vcf automatically if available) # colnames(vo) = i m0 = genCountMatrixFromVcf(mygenome, vo) mut = rbind(mut, m0) } dim(mut) # matrix with all samples
The opportunity matrix can also be generated from the reference genome (hg19 in the following case):
library(rtracklayer) target_regions <- import(con="/path/to/a/target.bed", format="bed") opp <- genOpportunityFromGenome(BSgenome.Hsapiens.UCSC.hg19, target_regions, nsamples=nrow(mut))
Where target.bed is a bed file
containing the genomic regions analyzed by the variant caller.
If a BSgenome is not available for your genome, you can use a fasta file:
library(Rsamtools) # make sure /path/to/genome.fasta.fai exists ! # you can use "samtools faidx" command to create it mygenome <- FaFile("/path/to/genome.fasta") mut <- genCountMatrixFromVcf(mygenome, vcfobj) opp <- genOpportunityFromGenome(mygenome, target_regions)
By convention, the input file should be tab-delimited with sample names as row names and features as column names. Features should be refered in the format "base change:triplet", e.g. "C>A:TCG", as can be seen in the example below. Similarly, the opportunity matrix can be provided in a tab-delimited file with the same structure as the mutation counts file. An example of the required matrix format can be seen here.
This tutorial uses as input the 21 breast cancer dataset described in
Nik-Zainal et al 2012. For the sake of convenience this dataset is
included with the package and can be accessed by using the
system.file function:
mut <- read.table(system.file("extdata","21_breast_cancers.mutations.txt", package="signeR"), header=TRUE, check.names=FALSE) opp <- read.table(system.file("extdata","21_breast_cancers.opportunity.txt", package="signeR"))
signeR takes a count matrix as its only required parameter, but the user can provide an opportunity matrix as well. The algorithm allows the assessment of the number of signatures by three options, as follows.
signatures <- signeR(M=mut, Opport=opp)
signatures <- signeR(M=mut, Opport=opp, nlim=c(2,11))
signatures <- signeR(M=mut, Opport=opp, nsig=5, main_eval=100, EM_eval=50, EMit_lim=20)
The parameters testing_burn and testing_eval
control the number of iterations used to estimate the number of signatures
(default value is 1000 for both parameters). There are other
arguments that may be passed on to signeR. Please have a look at signeR's
manual, issued by typing help(signeR).
Whenever signeR is left to decide which number of signatures is optimal, it will search for the rank Nsig that maximizes the median Bayesian Information Criterion (BIC). After the processing is done, this information can be plotted by the following command:
BICboxplot(signatures)
Boxplot of BIC values, showing that the optimal number of signatures for this dataset is 5.
The following instruction plots the MCMC sampled paths for each entry of the signature matrix P and their exposures, i.e. the E matrix. Only post-burnin paths are available for plotting. Those plots are useful for checking if entries have leveled off, reflecting the sampler convergence.
Paths(signatures$SignExposures)
Each plot shows the entries and exposures of one signature along sampler iterations.
Once the processing is done, the estimated signatures can be displayed in two charts, described below.
Signatures can be visualized in a barplot by issuing the following command:
SignPlot(signatures$SignExposures)
Signatures barplot with error bars reflecting the sample percentiles 0.05, 0.25, 0.75, and 0.95 for each entry.
Estimated signatures can also be visualized in a heatmap, generated by the following command:
SignHeat(signatures$SignExposures)
Heatmap showing the entries of each signature.
The relative contribution of each signature to the inspected genomes can be displayed in several ways. signeR currently provides three alternatives.
The levels of exposure to each signature in all genome samples can also be plotted:
ExposureBoxplot(signatures$SignExposures)
The relative contribution of the signatures to the mutations found on each genome sample can also be visualized by the following command:
ExposureBarplot(signatures$SignExposures)
Barplot showing the contributions of the signatures to genome samples mutation counts.
The levels of exposure to each signature can also be plotted in a heatmap:
ExposureHeat(signatures$SignExposures)
Heatmap showing the exposures for each genome sample. Samples are grouped according to their levels of exposure to the signatures, as can be seen in the dendrogram on the left.
signeR can highlight signatures that are differentially active among
previously defined groups of samples. To perform this task signeR needs
a vector of group labels. In this example the samples were divided according to
germline mutations at BRCA genes: groups wt,
BRCA1+ and BRCA2+, taken from the description of the
21 breast cancer data set.
# group labels, respective to each row of the mutation count matrix BRCA_labels <- c("wt","BRCA1+","BRCA2+","BRCA1+","BRCA2+","BRCA1+","BRCA1+", "wt","wt","wt","wt","BRCA1+","wt","BRCA2+","BRCA2+","wt","wt","wt", "wt","wt","wt") diff_exposure <- DiffExp(signatures$SignExposures, labels=BRCA_labels)
Top chart: DES plot showing that the BRCA+ samples were significantly more exposed to signatures S3, S4 and S5. Bottom chart: plots showing the significant differences found when groups are compared against each other. These last plots are generated only when there are more than two groups in the analysis and any signature is found to be differentially active. Groups marked by the same letter are not significantly different for the corresponding signature.
The Pvquant vector holds the pvalues quantile of the test for
each signature (by default, the 0.5 quantile, i.e. the median). The logarithms
of those are considered the Differential Exposure Scores (DES). Signatures with
Pvquant values below the cutoff, 0.05 by default, are considered as
differentially exposed.
# pvalues diff_exposure$Pvquant
## [1] 0.7452909713 0.4481633903 0.0017156464 0.0006192429 0.0303415234
The MostExposed vector contains the name of the group where each
differentially exposed signature showed the highest levels of activity.
# most exposed group diff_exposure$MostExposed
## [1] NA NA "BRCA2+" "BRCA1+" "BRCA1+"
signeR can also classify unlabeled samples based on the given labels.
In order to do this, those samples must correspond to NA values in
the labels vector and the Classify function can be used to assign them to one of
the defined groups. This example uses the sample labels defined in the DES
analysis performed previously.
# note that BRCA_labels [15],[20] and [21] are set to NA BRCA_labels <- c("wt","BRCA+","BRCA+","BRCA+","BRCA+","BRCA+","BRCA+","wt","wt", "wt","wt","BRCA+","wt","BRCA+",NA,"wt","wt","wt","wt",NA,NA) Class <- Classify(signatures$SignExposures, labels=BRCA_labels)
Barplot showing the relative frequencies of assignment of each unlabeled sample to the selected group.
# Final assignments Class$class
## [1] "BRCA+" "wt" "wt"
# Relative frequencies of assignment to selected groups Class$freq
## PD4116a PD4199a PD4248a ## 1 1 1
# All assigment frequencies Class$allfreqs
## PD4116a-BRCA+ PD4199a-wt PD4248a-wt ## BRCA+ 100 0 0 ## wt 0 100 100
citation("signeR")
##
## Rafael A. Rosales, Rodrigo D. Drummond, Renan Valieris, Emmanuel
## Dias-Neto, and Israel T. da Silva (2016): signeR: An empirical
## Bayesian approach to mutational signature discovery.
## Bioinformatics September 1, 2016
## doi:10.1093/bioinformatics/btw572
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {signeR: An empirical Bayesian approach to mutational signature discovery.},
## author = {Rafael A. Rosales and Rodrigo D. Drummond and Renan Valieris and Emmanuel Dias-Neto and Israel T. da Silva},
## year = {2016},
## journal = {Bioinformatics},
## doi = {10.1093/bioinformatics/btw572},
## }
OS X users might get compilation errors similar to these:
ld: warning: directory not found for option '-L/usr/local/lib/x86_64'
ld: library not found for -lgfortran
ld: library not found for -lquadmath
This problem arises when the machine is missing gfortran libraries necessary to compile RcppArmadillo and signeR. To install the missing libraries, execute these lines on a terminal:
curl -O http://r.research.att.com/libs/gfortran-4.8.2-darwin13.tar.bz2
sudo tar fvxz gfortran-4.8.2-darwin13.tar.bz2 -C /
For more information see this post and the Rcpp FAQ, section 2.16.
Some packages that signeR depends on requires that 3rd party library headers be installed. If you see errors like:
Error: checking for curl-config... no
Cannot find curl-configIt means you need to install these headers with your package manager. For example on ubuntu:
$ sudo apt-get install libcurl4-openssl-devL. B. Alexandrov, S. Nik-Zainal, D. C. Wedge, P. J. Campbell, and M. R. Stratton. Deciphering Signatures of Mutational Processes Operative in Human Cancer. Cell Reports, 3(1):246-259, Jan. 2013. [ DOI ]
A. Fischer, C. J. Illingworth, P. J. Campbell, and V. Mustonen. EMu: probabilistic inference of mutational processes and their localization in the cancer genome. Genome biology, 14(4):R39, Apr. 2013. [ DOI ]
sessionInfo()
## R version 3.6.0 (2019-04-26) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 18.04.2 LTS ## ## Matrix products: default ## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so ## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 parallel stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] signeR_1.10.0 NMF_0.21.0 ## [3] bigmemory_4.5.33 cluster_2.0.9 ## [5] rngtools_1.3.1.1 pkgmaker_0.27 ## [7] registry_0.5-1 VariantAnnotation_1.30.0 ## [9] Rsamtools_2.0.0 Biostrings_2.52.0 ## [11] XVector_0.24.0 SummarizedExperiment_1.14.0 ## [13] DelayedArray_0.10.0 BiocParallel_1.18.0 ## [15] matrixStats_0.54.0 Biobase_2.44.0 ## [17] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0 ## [19] IRanges_2.18.0 S4Vectors_0.22.0 ## [21] BiocGenerics_0.30.0 knitr_1.22 ## ## loaded via a namespace (and not attached): ## [1] httr_1.4.0 bit64_0.9-7 ## [3] foreach_1.4.4 assertthat_0.2.1 ## [5] highr_0.8 blob_1.1.1 ## [7] BSgenome_1.52.0 GenomeInfoDbData_1.2.1 ## [9] progress_1.2.0 pillar_1.3.1 ## [11] RSQLite_2.1.1 lattice_0.20-38 ## [13] glue_1.3.1 digest_0.6.18 ## [15] RColorBrewer_1.1-2 colorspace_1.4-1 ## [17] Matrix_1.2-17 plyr_1.8.4 ## [19] XML_3.98-1.19 pkgconfig_2.0.2 ## [21] bibtex_0.4.2 biomaRt_2.40.0 ## [23] zlibbioc_1.30.0 purrr_0.3.2 ## [25] xtable_1.8-4 scales_1.0.0 ## [27] tibble_2.1.1 ggplot2_3.1.1 ## [29] withr_2.1.2 GenomicFeatures_1.36.0 ## [31] lazyeval_0.2.2 mime_0.6 ## [33] magrittr_1.5 crayon_1.3.4 ## [35] memoise_1.1.0 evaluate_0.13 ## [37] doParallel_1.0.14 class_7.3-15 ## [39] PMCMR_4.3 tools_3.6.0 ## [41] prettyunits_1.0.2 hms_0.4.2 ## [43] gridBase_0.4-7 stringr_1.4.0 ## [45] munsell_0.5.0 AnnotationDbi_1.46.0 ## [47] compiler_3.6.0 rlang_0.3.4 ## [49] nloptr_1.2.1 grid_3.6.0 ## [51] RCurl_1.95-4.12 iterators_1.0.10 ## [53] bigmemory.sri_0.1.3 bitops_1.0-6 ## [55] gtable_0.3.0 codetools_0.2-16 ## [57] DBI_1.0.0 markdown_0.9 ## [59] reshape2_1.4.3 R6_2.4.0 ## [61] GenomicAlignments_1.20.0 dplyr_0.8.0.1 ## [63] rtracklayer_1.44.0 bit_1.1-14 ## [65] stringi_1.4.3 Rcpp_1.0.1 ## [67] tidyselect_0.2.5 xfun_0.6
print(names(dev.cur()))
## [1] "png"