selectKSigs: a package for selecting the number of mutational signatures

In selectKSigs: Selecting the number of mutational signatures using a perplexity-based measure and cross-validation

library(BiocStyle)

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.

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.

Installing and loading the package {#installation}

Installation

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 r CRANpkg("BiocManager").]

You can also download the newest version from the GitHub using devtools:

devtools::install_github("USCbiostats/selectKSigs")

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.

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).

Workflow

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)
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"))

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.

set.seed(5)
results <- cv_PMSignature(G, Kfold = 3, nRep = 3, Klimit = 7)
print(results)

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")

Session info

Here is the output of sessionInfo() for reproducibility in the future.

sessionInfo()

selectKSigs documentation built on Nov. 8, 2020, 5:38 p.m.