## ----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('

',options$fig.cap,"

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