## ----results='hide', echo=FALSE, message=FALSE, warning=FALSE-------
set.seed(1)

options(width = 70)

library(knitr)

style_sheet = "bioc.css"
style = if(file.exists(style_sheet)) {
    paste(readLines(style_sheet), collapse = "\n")
}

opts_knit$set(self.contained = TRUE,
              upload.fun = image_uri,
              header = c(highlight = style))

opts_chunk$set(comment = "  ",
               fig.path = "",
               fig.align = "center",
               out.width = "50%",
               dpi = 300,
               indent = 10,
               cache = FALSE,
               cache.path = "../cache")

knit_hooks$set(fig.cap = function(before, options, envir) {
    if(!before) {
        paste('<p class="caption">',options$fig.cap,"</p>",sep="")
    }
})

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

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

## ----load_data_package, results='hide',message=FALSE----------------
library(SomaticCancerAlterations)
library(BSgenome.Hsapiens.UCSC.hg19)

## ----sca_metadata---------------------------------------------------
sca_metadata = scaMetadata()

print(sca_metadata)

## ----sca_to_vranges-------------------------------------------------
sca_vr = scaSNVRanges()

head(sca_vr, 3)

## ----sca_study_table------------------------------------------------
sort(table(sca_vr$study), decreasing = TRUE)

## ----sca_vr_to_motifs-----------------------------------------------
sca_motifs = mutationContext(sca_vr, BSgenome.Hsapiens.UCSC.hg19, unify = TRUE)

## ----sca_motif_occurrence-------------------------------------------
sca_mm = motifMatrix(sca_motifs, group = "study", normalize = TRUE)

head(round(sca_mm, 4))

## ----sca_mutation_spectrum, fig.cap='Mutation spectrum over studies'----
plotMutationSpectrum(sca_motifs, "study")

## ----sca_nmf_pca----------------------------------------------------
n_sigs = 5

sigs_nmf = identifySignatures(sca_mm, n_sigs, nmfDecomposition)

sigs_pca = identifySignatures(sca_mm, n_sigs, pcaDecomposition)

## ----sca_explore_nmf------------------------------------------------
sigs_nmf

## ----sca_explore_pca------------------------------------------------
sigs_pca

## ----sca_plot_nmf_signatures_map, fig.cap='Composition of somatic signatures estimated with the NMF, represented as a heatmap.'----
plotSignatureMap(sigs_nmf) + ggtitle("Somatic Signatures: NMF - Heatmap")

## ----sca_plot_nmf_signatures, fig.cap='Composition of somatic signatures estimated with the NMF, represented as a barchart.'----
plotSignatures(sigs_nmf) + ggtitle("Somatic Signatures: NMF - Barchart")

## -------------------------------------------------------------------
plotObservedSpectrum(sigs_nmf)

## -------------------------------------------------------------------
plotFittedSpectrum(sigs_nmf)

## ----sca_plot_nmf_samples_map, fig.cap='Occurrence of signatures estimated with the NMF, represented as a heatmap.'----
plotSampleMap(sigs_nmf)

## ----sca_plot_nmf_samples, fig.cap='Occurrence of signatures estimated with the NMF, represented as a barchart.'----
plotSamples(sigs_nmf)

## ----sca_plot_pca_signatures_map, fig.cap='Composition of somatic signatures estimated with the PCA, represented as a heatmap.'----
plotSignatureMap(sigs_pca) + ggtitle("Somatic Signatures: PCA - Heatmap")

## ----sca_plot_pca_signatures, fig.cap='Composition of somatic signatures estimated with the PCA, represented as a barchart.'----
plotSignatures(sigs_pca) + ggtitle("Somatic Signatures: PCA - Barchart")

## -------------------------------------------------------------------
plotObservedSpectrum(sigs_pca)

## -------------------------------------------------------------------
plotFittedSpectrum(sigs_pca)

## ----sva_batch_not_run, eval=FALSE----------------------------------
#  library(sva)
#  
#  df = as(sca_metadata, "data.frame") ## sample x covariable
#  pheno = data.frame(s = unlist(df[ ,"Sequence_Source"]), c = unlist(df[ ,"Cancer_Type"]))
#  rownames(pheno) = gsub("(.*)_.*", "\\1", rownames(pheno))
#  mod = model.matrix(~ s + c, data = pheno)
#  mod0 = model.matrix(~ c, data = pheno)
#  
#  sv = sva(sca_mm, mod, mod0, method = "irw")

## ----kmer_chr1, eval=FALSE------------------------------------------
#  k = 3
#  n = 1e5
#  chrs = "chr1"
#  
#  chr1_ranges = as(seqinfo(BSgenome.Hsapiens.UCSC.hg19), "GRanges")
#  chr1_ranges = keepSeqlevels(chr1_ranges, chrs)
#  
#  k3_chr1 = kmerFrequency(BSgenome.Hsapiens.UCSC.hg19, n, k, chr1_ranges)
#  
#  k3_chr1

## ----normalize_motifs, eval=FALSE-----------------------------------
#  head(sca_mm)
#  
#  data(kmers)
#  norms = k3wg / k3we
#  head(norms)
#  
#  sca_norm = normalizeMotifs(sca_mm, norms)
#  
#  head(sca_norm)

## -------------------------------------------------------------------
citation("SomaticSignatures")

## ----eval=FALSE-----------------------------------------------------
#  source("http://bioconductor.org/biocLite.R")
#  biocLite("SomaticSignatures")

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

## -------------------------------------------------------------------
vr = VRanges(
    seqnames = "chr1",
    ranges = IRanges(start = 1000, width = 1),
    ref = "A",
    alt = "C")

vr

## ----echo=FALSE, results='markup'-----------------------------------
sessionInfo()