library(BiocStyle)
The R package HiLDA is developed under the Bayesian framework to select the number of mutational signatures based on perplexity and cross-validation. The mutation signature is defined based on the independent model proposed by Shiraishi's et al.
Shiraishi et al. A simple model-based approach to inferring and visualizing cancer mutation signatures, bioRxiv, doi: http://dx.doi.org/10.1101/019901.
Zhi Yang, Paul Marjoram, Kimberly D. Siegmund. Selecting the number of mutational signatures using a perplexity-based measure and cross-validation.
selectKSigs requires several CRAN and Bioconductor R packages to be installed. Dependencies are usually handled automatically, when installing the package using the following commands:
install.packages("BiocManager") BiocManager::install("selectKSigs")
[NOTE: Ignore the first line if you already have installed the
r CRANpkg("BiocManager")
.]
You can also download the newest version from the GitHub using devtools:
devtools::install_github("USCbiostats/selectKSigs")
selectKSigs
is a package built on some basic functions from pmsignature
including how to read the input data. Here is an example from pmsignature
on
the input data, mutation features are elements used for categorizing mutations
such as:
sample1 chr1 100 A C sample1 chr1 200 A T sample1 chr2 100 G T sample2 chr1 300 T C sample3 chr3 400 T C
Here, inputFile is the path for the input file. numBases is the number of flanking bases to consider including the central base (if you want to consider two 5' and 3' bases, then set 5). Also, you can add transcription direction information using trDir. numSig sets the number of mutation signatures estimated from the input data. You will see a warning message on some mutations are being removed.
library(HiLDA) library(tidyr) library(ggplot2) library(dplyr) inputFile <- system.file("extdata/esophageal.mp.txt.gz", package="HiLDA") G <- hildaReadMPFile(inputFile, numBases=5, trDir=TRUE)
Also, we also provided a small simulated dataset which contains 10 mutational catalogs and used it for demonstrating the key functions in selectKSigs. We start with loading the sample dataset G stored as extdata/sample.rdata.
library(selectKSigs) load(system.file("extdata/sample.rdata", package = "selectKSigs"))
After we read in the sample data G, we can run the process from selectKSigs. Here, we specify the inputG as G, the number of cross-validation folds, kfold to be 3, the number of replications, nRep, to be 3, and the upper limit of the K values for exploration to be 7.
set.seed(5) results <- cv_PMSignature(G, Kfold = 3, nRep = 3, Klimit = 7) print(results)
After we obtained the results, we can plot each measure by the range of K values that were refitted during the calculation. The optimal value of K is achieved at its minimum value highlighted in grey.
results$Kvalue <- seq_len(nrow(results)) + 1 results_df <- gather(results, Method, value, -Kvalue) %>% group_by(Method) %>% mutate(xmin = which.min(value) + 1 - 0.1, xmax = which.min(value) + 1 + 0.1) ggplot(results_df) + geom_point(aes(x = Kvalue, y = value, color = Method), size = 2) + facet_wrap(~ Method, scales = "free") + geom_rect(mapping = aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf), fill = 'grey', alpha = 0.05) + theme_bw()+ xlab("Number of signatures")
Here is the output of sessionInfo()
for reproducibility in the future.
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.