1 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: ggplot2
library(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, union
inputFile <- 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.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.11-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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] selectKSigs_1.0.0 dplyr_0.8.5       tidyr_1.0.2       HiLDA_1.2.0      
## [5] ggplot2_3.3.0     BiocStyle_2.16.0 
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6                           
##  [2] matrixStats_0.56.0                     
##  [3] bit64_0.9-7                            
##  [4] progress_1.2.2                         
##  [5] httr_1.4.1                             
##  [6] GenomeInfoDb_1.24.0                    
##  [7] tools_4.0.0                            
##  [8] R6_2.4.1                               
##  [9] R2WinBUGS_2.1-21                       
## [10] DBI_1.1.0                              
## [11] BiocGenerics_0.34.0                    
## [12] colorspace_1.4-1                       
## [13] withr_2.2.0                            
## [14] tidyselect_1.0.0                       
## [15] prettyunits_1.1.1                      
## [16] bit_1.1-15.2                           
## [17] curl_4.3                               
## [18] compiler_4.0.0                         
## [19] Biobase_2.48.0                         
## [20] DelayedArray_0.14.0                    
## [21] rtracklayer_1.48.0                     
## [22] labeling_0.3                           
## [23] bookdown_0.18                          
## [24] scales_1.1.0                           
## [25] askpass_1.1                            
## [26] rappdirs_0.3.1                         
## [27] stringr_1.4.0                          
## [28] digest_0.6.25                          
## [29] Rsamtools_2.4.0                        
## [30] rmarkdown_2.1                          
## [31] BSgenome.Hsapiens.UCSC.hg19_1.4.3      
## [32] XVector_0.28.0                         
## [33] pkgconfig_2.0.3                        
## [34] htmltools_0.4.0                        
## [35] dbplyr_1.4.3                           
## [36] BSgenome_1.56.0                        
## [37] rlang_0.4.5                            
## [38] RSQLite_2.2.0                          
## [39] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [40] farver_2.0.3                           
## [41] BiocParallel_1.22.0                    
## [42] gtools_3.8.2                           
## [43] RCurl_1.98-1.2                         
## [44] magrittr_1.5                           
## [45] GenomeInfoDbData_1.2.3                 
## [46] Matrix_1.2-18                          
## [47] Rcpp_1.0.4.6                           
## [48] munsell_0.5.0                          
## [49] S4Vectors_0.26.0                       
## [50] abind_1.4-5                            
## [51] lifecycle_0.2.0                        
## [52] stringi_1.4.6                          
## [53] yaml_2.2.1                             
## [54] SummarizedExperiment_1.18.0            
## [55] zlibbioc_1.34.0                        
## [56] BiocFileCache_1.12.0                   
## [57] grid_4.0.0                             
## [58] blob_1.2.1                             
## [59] parallel_4.0.0                         
## [60] forcats_0.5.0                          
## [61] crayon_1.3.4                           
## [62] lattice_0.20-41                        
## [63] Biostrings_2.56.0                      
## [64] cowplot_1.0.0                          
## [65] GenomicFeatures_1.40.0                 
## [66] hms_0.5.3                              
## [67] magick_2.3                             
## [68] knitr_1.28                             
## [69] pillar_1.4.3                           
## [70] GenomicRanges_1.40.0                   
## [71] boot_1.3-25                            
## [72] biomaRt_2.44.0                         
## [73] stats4_4.0.0                           
## [74] XML_3.99-0.3                           
## [75] R2jags_0.6-1                           
## [76] glue_1.4.0                             
## [77] evaluate_0.14                          
## [78] BiocManager_1.30.10                    
## [79] vctrs_0.2.4                            
## [80] gtable_0.3.0                           
## [81] openssl_1.4.1                          
## [82] purrr_0.3.4                            
## [83] assertthat_0.2.1                       
## [84] xfun_0.13                              
## [85] coda_0.19-3                            
## [86] rjags_4-10                             
## [87] tibble_3.0.1                           
## [88] GenomicAlignments_1.24.0               
## [89] AnnotationDbi_1.50.0                   
## [90] memoise_1.1.0                          
## [91] IRanges_2.22.0                         
## [92] ellipsis_0.3.0