inst/doc/signeR-vignette.R

In signeR: Empirical Bayesian approach to mutational signature discovery

## ----results='hide', echo=FALSE, message=FALSE, warning=FALSE-----------------
set.seed(42)
library(knitr)
opts_knit$set(
    self.contained = TRUE,
    upload.fun = image_uri
)
opts_chunk$set(
    dev = 'png',
    dpi = 300,
    out.width = "700px",
    out.height = "700px"
)

## ----eval=FALSE---------------------------------------------------------------
#  install.packages("BiocManager")
#  BiocManager::install("signeR")

## ----results='hide', message=FALSE--------------------------------------------
library(signeR)

## ----eval=FALSE---------------------------------------------------------------
#  library(VariantAnnotation)
#  
#  # BSgenome, equivalent to the one used on the variant call
#  library(BSgenome.Hsapiens.UCSC.hg19)
#  
#  vcfobj <- readVcf("/path/to/a/file.vcf", "hg19")
#  mut <- genCountMatrixFromVcf(BSgenome.Hsapiens.UCSC.hg19, vcfobj)

## ----eval=FALSE---------------------------------------------------------------
#  mut = matrix(ncol=96,nrow=0)
#  for(i in vcf_files) {
#      vo = readVcf(i, "hg19")
#      # sample name (should pick up from the vcf automatically if available)
#      # colnames(vo) = i
#      m0 = genCountMatrixFromVcf(mygenome, vo)
#      mut = rbind(mut, m0)
#  }
#  dim(mut) # matrix with all samples

## ----eval=FALSE---------------------------------------------------------------
#  library(rtracklayer)
#  
#  target_regions <- import(con="/path/to/a/target.bed", format="bed")
#  opp <- genOpportunityFromGenome(BSgenome.Hsapiens.UCSC.hg19,
#      target_regions, nsamples=nrow(mut))

## ----eval=FALSE---------------------------------------------------------------
#  library(Rsamtools)
#  
#  # make sure /path/to/genome.fasta.fai exists !
#  # you can use "samtools faidx" command to create it
#  mygenome <- FaFile("/path/to/genome.fasta")
#  
#  mut <- genCountMatrixFromVcf(mygenome, vcfobj)
#  opp <- genOpportunityFromGenome(mygenome, target_regions)
#  

## -----------------------------------------------------------------------------
mut <- read.table(system.file("extdata","21_breast_cancers.mutations.txt",
    package="signeR"), header=TRUE, check.names=FALSE)
opp <- read.table(system.file("extdata","21_breast_cancers.opportunity.txt",
    package="signeR"))

## ----eval=FALSE---------------------------------------------------------------
#  signatures <- signeR(M=mut, Opport=opp)

## ----eval=FALSE---------------------------------------------------------------
#  signatures <- signeR(M=mut, Opport=opp, nlim=c(2,11))

## ----results='hide', message=FALSE--------------------------------------------
signatures <- signeR(M=mut, Opport=opp, nsig=5, main_eval=100, EM_eval=50, EMit_lim=20)

## ----eval=FALSE---------------------------------------------------------------
#  BICboxplot(signatures)

## ----echo=FALSE, results='asis'-----------------------------------------------
cat(sprintf("<img src=\"%s\" />\n",image_uri("Model_selection_BICs.png")))

## -----------------------------------------------------------------------------
Paths(signatures$SignExposures)

## -----------------------------------------------------------------------------
SignPlot(signatures$SignExposures)

## -----------------------------------------------------------------------------
SignHeat(signatures$SignExposures)

## -----------------------------------------------------------------------------
ExposureBoxplot(signatures$SignExposures)

## -----------------------------------------------------------------------------
ExposureBarplot(signatures$SignExposures)

## -----------------------------------------------------------------------------
ExposureHeat(signatures$SignExposures)

## -----------------------------------------------------------------------------
# group labels, respective to each row of the mutation count matrix
BRCA_labels <- c("wt","BRCA1+","BRCA2+","BRCA1+","BRCA2+","BRCA1+","BRCA1+",
    "wt","wt","wt","wt","BRCA1+","wt","BRCA2+","BRCA2+","wt","wt","wt",
    "wt","wt","wt")

diff_exposure <- DiffExp(signatures$SignExposures, labels=BRCA_labels)

## -----------------------------------------------------------------------------
# pvalues
diff_exposure$Pvquant

## -----------------------------------------------------------------------------
# most exposed group
diff_exposure$MostExposed

## -----------------------------------------------------------------------------
# note that BRCA_labels [15],[20] and [21] are set to NA
BRCA_labels <- c("wt","BRCA+","BRCA+","BRCA+","BRCA+","BRCA+","BRCA+","wt","wt",
    "wt","wt","BRCA+","wt","BRCA+",NA,"wt","wt","wt","wt",NA,NA)

Class <- Classify(signatures$SignExposures, labels=BRCA_labels)

## -----------------------------------------------------------------------------
# Final assignments
Class$class

# Relative frequencies of assignment to selected groups
Class$freq

# All assigment frequencies
Class$allfreqs

## -----------------------------------------------------------------------------
citation("signeR")

## -----------------------------------------------------------------------------
sessionInfo()

## -----------------------------------------------------------------------------
print(names(dev.cur()))