selectKSigs: a package for selecting the number of mutational signatures
2021-05-19
Abstract
Instructions on using selectKSigs on selecting the number of mutational signatures using a perplexity-based measure and cross-validation1 Introduction
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.Â
1.1 Paper
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.
2 Installing and loading the package
2.1 Installation
2.1.1 Bioconductor
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 BiocManager.]
You can also download the newest version from the GitHub using devtools:
devtools::install_github("USCbiostats/selectKSigs")3 Input data
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:
- 6 substitutions (C>A, C>G, C>T, T>A, T>C and T>G)
- 2 flanking bases (A, C, G and T)
- transcription direction.
3.1 Mutation Position Format
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 - The 1st column shows the name of samples
- The 2nd column shows the name of chromosome
- The 3rd column shows the coordinate in the chromosome
- The 4th column shows the reference base (A, C, G, or T).
- The 5th colum shows the alternate base (A, C, G, or T).
4 Workflow
4.1 Get input data
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)## Loading required package: ggplot2library(tidyr)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioninputFile <- system.file("extdata/esophageal.mp.txt.gz", package="HiLDA")
G <- hildaReadMPFile(inputFile, numBases=5, trDir=TRUE)## Warning in hildaReadMPFile(inputFile, numBases = 5, trDir = TRUE): Out of
## 24861 mutations, we could obtaintranscription direction information for 24728
## mutation. Other mutations are removed.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"))4.1.1 Perform the selecting process
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.
4.1.2 Visualizing the 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")
4.2 Session info
Here is the output of sessionInfo() for reproducibility in the future.
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.13-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.13-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] selectKSigs_1.4.0 dplyr_1.0.6 tidyr_1.1.3 HiLDA_1.6.0
## [5] ggplot2_3.3.3 BiocStyle_2.20.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7
## [2] matrixStats_0.58.0
## [3] bit64_4.0.5
## [4] filelock_1.0.2
## [5] progress_1.2.2
## [6] httr_1.4.2
## [7] GenomeInfoDb_1.28.0
## [8] tools_4.1.0
## [9] bslib_0.2.5.1
## [10] utf8_1.2.1
## [11] R6_2.5.0
## [12] R2WinBUGS_2.1-21
## [13] DBI_1.1.1
## [14] BiocGenerics_0.38.0
## [15] colorspace_2.0-1
## [16] withr_2.4.2
## [17] prettyunits_1.1.1
## [18] tidyselect_1.1.1
## [19] curl_4.3.1
## [20] bit_4.0.4
## [21] compiler_4.1.0
## [22] Biobase_2.52.0
## [23] DelayedArray_0.18.0
## [24] labeling_0.4.2
## [25] rtracklayer_1.52.0
## [26] bookdown_0.22
## [27] sass_0.4.0
## [28] scales_1.1.1
## [29] rappdirs_0.3.3
## [30] stringr_1.4.0
## [31] digest_0.6.27
## [32] Rsamtools_2.8.0
## [33] rmarkdown_2.8
## [34] BSgenome.Hsapiens.UCSC.hg19_1.4.3
## [35] XVector_0.32.0
## [36] pkgconfig_2.0.3
## [37] htmltools_0.5.1.1
## [38] MatrixGenerics_1.4.0
## [39] highr_0.9
## [40] dbplyr_2.1.1
## [41] fastmap_1.1.0
## [42] BSgenome_1.60.0
## [43] rlang_0.4.11
## [44] RSQLite_2.2.7
## [45] farver_2.1.0
## [46] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [47] jquerylib_0.1.4
## [48] BiocIO_1.2.0
## [49] generics_0.1.0
## [50] jsonlite_1.7.2
## [51] gtools_3.8.2
## [52] BiocParallel_1.26.0
## [53] RCurl_1.98-1.3
## [54] magrittr_2.0.1
## [55] GenomeInfoDbData_1.2.6
## [56] Matrix_1.3-3
## [57] Rcpp_1.0.6
## [58] munsell_0.5.0
## [59] S4Vectors_0.30.0
## [60] fansi_0.4.2
## [61] abind_1.4-5
## [62] lifecycle_1.0.0
## [63] stringi_1.6.2
## [64] yaml_2.2.1
## [65] SummarizedExperiment_1.22.0
## [66] zlibbioc_1.38.0
## [67] BiocFileCache_2.0.0
## [68] grid_4.1.0
## [69] blob_1.2.1
## [70] parallel_4.1.0
## [71] forcats_0.5.1
## [72] crayon_1.4.1
## [73] lattice_0.20-44
## [74] Biostrings_2.60.0
## [75] cowplot_1.1.1
## [76] GenomicFeatures_1.44.0
## [77] hms_1.1.0
## [78] KEGGREST_1.32.0
## [79] magick_2.7.2
## [80] knitr_1.33
## [81] pillar_1.6.1
## [82] GenomicRanges_1.44.0
## [83] boot_1.3-28
## [84] rjson_0.2.20
## [85] biomaRt_2.48.0
## [86] stats4_4.1.0
## [87] XML_3.99-0.6
## [88] R2jags_0.6-1
## [89] glue_1.4.2
## [90] evaluate_0.14
## [91] BiocManager_1.30.15
## [92] vctrs_0.3.8
## [93] png_0.1-7
## [94] gtable_0.3.0
## [95] purrr_0.3.4
## [96] assertthat_0.2.1
## [97] cachem_1.0.5
## [98] xfun_0.23
## [99] restfulr_0.0.13
## [100] coda_0.19-4
## [101] rjags_4-10
## [102] tibble_3.1.2
## [103] GenomicAlignments_1.28.0
## [104] AnnotationDbi_1.54.0
## [105] memoise_2.0.0
## [106] IRanges_2.26.0
## [107] ellipsis_0.3.2