---
title: "selectKSigs: a package for selecting the number of mutational signatures"
author: 
- name: Zhi Yang
  affiliation: Department of Preventive Medicine, 
    University of Southern California, Los Angeles, USA 
  email: zyang895@gmail.com
date: "`r Sys.Date()`"
vignette: |
  %\VignetteIndexEntry{An introduction to HiLDA}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
output:
  BiocStyle::html_document:
    toc_float: true
  BiocStyle::pdf_document: default
abstract: | 
  Instructions on using _selectKSigs_ on selecting the number of mutational 
  signatures using a perplexity-based measure and cross-validation

---

```{r style, echo = FALSE, results = 'asis'}
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](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. 

```{r}
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.
```{r}
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. 
```{r include=FALSE}
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. 

```{r}
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. 

```{r}
sessionInfo()
```