inst/doc/Rariant-vignette.R

In Rariant: Identification and Assessment of Single Nucleotide Variants through Shifts in Non-Consensus Base Call Frequencies

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

options(width = 70)

style_file = "bioc.css"

library(knitr)

style = if(file.exists(style_file))
            paste(readLines(style_file), 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 = "100%",
               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="")
    }
})

## ----ci_cases_plot, echo=FALSE, message=FALSE, fig.width=14, fig.height=7, fig.cap='Illustrative cases of confidence intervals for somatic variant frequency estimates'----
library(ggplot2)

df = data.frame(
    x = factor(rep(c(""), times = 6)),
    case = factor(1:6),
    d = c(0.5, 0.3, -0.525, 0, 0.2, 0),
    cil = c(0.45, 0.2, -0.60, -0.05, -0.3, -1),
    ciu = c(0.55, 0.4, -0.45, 0.05, 0.7, 1)
    )

p = ggplot(df) + geom_hline(aes(yintercept = 0), color = "darkgray") +
geom_pointrange(aes(x = x, y = d, ymin = cil, ymax = ciu), size = 1, color =
"black") + facet_grid( ~ case) + ylim(-1, 1) + theme_bw() +
theme(legend.position = "none") + xlab("") + ylab("pT - pC")

print(p)

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

library(GenomicRanges)
library(ggbio)

## -------------------------------------------------------------------
tp53_region = GRanges("chr17", IRanges(7565097, 7590856))

## -------------------------------------------------------------------
control_bam = system.file("extdata", "platinum", "control.bam", package = "Rariant", mustWork = TRUE)
test1_bam = system.file("extdata", "platinum", "test.bam", package = "Rariant", mustWork = TRUE)
test2_bam = system.file("extdata", "platinum", "test2.bam", package = "Rariant", mustWork = TRUE)
mix_bam = system.file("extdata", "platinum", "mix.bam", package = "Rariant", mustWork = TRUE)

## -------------------------------------------------------------------
v_test1 = rariant(test1_bam, control_bam, tp53_region, select = FALSE)
v_test2 = rariant(test2_bam, control_bam, tp53_region, select = FALSE)
v_mix = rariant(mix_bam, control_bam, tp53_region, select = FALSE)

## -------------------------------------------------------------------
v_all = GenomicRangesList(T1 = v_test1, T2 = v_test2, M = v_mix)

v_all = endoapply(v_all, updateCalls)

## ----fig.width=14, fig.height=7, fig.cap='Confidence intervals for simulation study'----
t_ci = tracks(lapply(v_all, plotConfidenceIntervals, color = "verdict")) + verdictColorScale()

print(t_ci)

## ----fig.width=14, fig.height=7, fig.cap='Confidence intervals for simulation study'----
t_ci = tracks(lapply(v_all, plotConfidenceIntervals, color = "eventType"))

print(t_ci)

## ----fig.width=14, fig.height=7, fig.cap='Abundance shifts for simulation study'----
t_rates = tracks(lapply(v_all, plotAbundanceShift))

print(t_rates)

## -------------------------------------------------------------------
z = filterCalls(v_all, verdict %in% c("present", "inbetween", "dontknow"))

elementNROWS(z)

## ----fig.width=7, fig.height=7--------------------------------------
evidenceHeatmap(z, fill = "d", color = "verdict") + verdictColorScale()

## ----echo=FALSE, message=FALSE, fig.width=14, fig.height=7, fig.cap='Illustrative cases of confidence intervals for somatic variant frequency estimates for two strands'----
library(ggplot2)

df = data.frame(
    x = factor(rep(c("A", "B"), times = 7)),
    case = factor(rep(c(5, 6, 7, 4, 1, 2, 3), each = 2)),
    dx = c(0.65, -0.65,  0.65,  0.20,  0.65,  0,  0.65,  0.55, 0.65, 0, 0.05, -0.05, -0.05, 0.05),
    cil = c(0.5, -0.8, 0.5, 0.1, 0.5, -0.2, 0.5, 0.4, 0.5, -0.7, -0.1, -0.2, -0.9, -0.8),
    ciu = c(0.8, -0.5, 0.8, 0.3, 0.8, 0.2, 0.8, 0.7, 0.8, 0.7, 0.2, 0.1, 0.8, 0.9),
    group = factor(c(rep("n", 2*3), rep("o", 2*4)))
    )

p = ggplot(df) + geom_hline(aes(yintercept = 0), color = "darkgray") + geom_hline(aes(yintercept = 0.6), color = "darkred", linetype = "dashed") + geom_pointrange(aes(x = x, y = dx, ymin = cil, ymax = ciu), size = 1, color = "black") + facet_wrap(~ case, nrow = 2) + ylim(-1, 1) +  theme_bw() + theme(legend.position = "none") + xlab("Strand") + ylab("Shift in non-consensus rate")

print(p)

## -------------------------------------------------------------------
## WGS
n1 = 30
n2 = 30
k1 = 0:(n1-1)
k2 = 0:(n2-1)
cl = 0.95
n_sample = 1e4

pars = expand.grid(k1 = k1, k2 = k2, n1 = n1, n2 = n2, conf_level = cl)

cp_ac = coverageProbability(pars, fun = acCi, n_sample = n_sample)

## ----fig.width=7, fig.height=7, out.width='60%', fig.cap='Coverage probabilities for whole-genome setting'----
p_ac = ggplot(cp_ac) + geom_tile(aes(x = k1, y = k2, fill = cp)) +
    scale_fill_gradient2(midpoint = 0.95, limits = c(0.9, 1)) +
    theme_bw() + xlab("kT") + ylab("kC")

print(p_ac)

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

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