## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", crop = NULL ## Related to https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016656.html ) ## ----vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE---------------- ## Track time spent on making the vignette startTime <- Sys.time() ## Bib setup library("RefManageR") ## Write bibliography information bib <- c( R = citation(), BiocStyle = citation("BiocStyle")[1], knitr = citation("knitr")[1], rmarkdown = citation("rmarkdown")[1], sessioninfo = citation("sessioninfo")[1], testthat = citation("testthat")[1], megadepth = citation("megadepth")[1], megadepthpaper = citation("megadepth")[2], bioconductor2015 = RefManageR::BibEntry( bibtype = "Article", key = "bioconductor2015", author = "Wolfgang Huber and Vincent J Carey and Robert Gentleman and Simon Anders and Marc Carlson and Benilton S Carvalho and Hector Corrada Bravo and Sean Davis and Laurent Gatto and Thomas Girke and Raphael Gottardo and Florian Hahne and Kasper D Hansen and Rafael A Irizarry and Michael Lawrence and Michael I Love and James MacDonald and Valerie Obenchain and Andrzej K Oleś and Hervé Pagès and Alejandro Reyes and Paul Shannon and Gordon K Smyth and Dan Tenenbaum and Levi Waldron and Martin Morgan", title = "Orchestrating high-throughput genomic analysis with Bioconductor", year = 2015, doi = "10.1038/nmeth.3252", journal = "Nature Methods", journaltitle = "Nat Methods" ), RefManageR = citation("RefManageR")[1], rtracklayer = citation("rtracklayer") ) ## ----"runtime", out.width="100%", echo = FALSE, fig.cap = "Timing results. Timing comparison for processing 1,000 genomic regions on a bigWig file that is available on the local disk or on a remote resource. We compared megadepth against rtracklayer and pyBigWig. megadepth is typically faster that these other software solutions for computing the mean coverage across a set of 1,000 input genomic regions. Check for more details."---- knitr::include_graphics("https://raw.githubusercontent.com/LieberInstitute/megadepth/master/analysis/md_rt_pybw_runtime.png") ## ----"install", eval = FALSE-------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # # BiocManager::install("megadepth") # # ## Check that you have a valid Bioconductor installation # BiocManager::valid() ## ----"citation"--------------------------------------------------------------- ## Citation info citation("megadepth") ## ----"start", message=FALSE--------------------------------------------------- library("megadepth") ## ----"install_software"------------------------------------------------------- ## Install the latest version of Megadepth install_megadepth(force = TRUE) ## ----"locate_example_bw"------------------------------------------------------ ## Next, we locate the example BigWig and annotation files example_bw <- system.file("tests", "test.bam.all.bw", package = "megadepth", mustWork = TRUE ) annotation_file <- system.file("tests", "testbw2.bed", package = "megadepth", mustWork = TRUE ) ## Where are they? example_bw annotation_file ## ----example------------------------------------------------------------------ ## We can then use megadepth to compute the coverage bw_cov <- get_coverage( example_bw, op = "mean", annotation = annotation_file ) bw_cov ## ----"interface_options"------------------------------------------------------ ## R-like interface ## that captures the standard output into R head(megadepth_shell(help = TRUE)) ## Command-like interface megadepth_cmd("--help") ## ----"show_help", echo = FALSE------------------------------------------------ x <- megadepth_shell(help = TRUE) cat(paste0(x, "\n")) ## ----"bam_to_bigwig", eval = !xfun::is_windows()------------------------------ ## Find the example BAM file example_bam <- system.file("tests", "test.bam", package = "megadepth", mustWork = TRUE ) ## Create the BigWig file ## Currently Megadepth does not support this on Windows example_bw <- bam_to_bigwig(example_bam, overwrite = TRUE) ## Path to the output files generated by bam_to_bigwig() example_bw ## ----"bam_to_bigwig_windows", eval = xfun::is_windows()----------------------- # ## On Windows, use the example bigWig file that is already included in # ## the R package # example_bw <- system.file("tests", "test.bam.all.bw", # package = "megadepth", mustWork = TRUE # ) ## ----"get_coverage"----------------------------------------------------------- ## Next, we locate the example annotation BED file annotation_file <- system.file("tests", "testbw2.bed", package = "megadepth", mustWork = TRUE ) ## Now we can compute the coverage bw_cov <- get_coverage(example_bw, op = "mean", annotation = annotation_file) bw_cov ## ----"rtracklayer_derfinder"-------------------------------------------------- ## Checking with derfinder and rtracklayer bed <- rtracklayer::import(annotation_file) ## The file needs a name names(example_bw) <- "example" ## Read in the base-pair coverage data regionCov <- derfinder::getRegionCoverage( regions = bed, files = example_bw, verbose = FALSE ) ## Summarize the base-pair coverage data. ## Note that we have to round the mean to make them comparable. testthat::expect_equivalent( round(sapply(regionCov[c(1, 3:4, 2)], function(x) mean(x$value)), 2), bw_cov$score, ) ## If we compute the sum, there's no need to round testthat::expect_equivalent( sapply(regionCov[c(1, 3:4, 2)], function(x) sum(x$value)), get_coverage(example_bw, op = "sum", annotation = annotation_file)$score, ) ## ----"bam_to_junctions"------------------------------------------------------- ## Find the example BAM file example_bam <- system.file("tests", "test.bam", package = "megadepth", mustWork = TRUE ) ## Run bam_to_junctions() example_jxs <- bam_to_junctions(example_bam, overwrite = TRUE) ## Path to the output file generated by bam_to_junctions() example_jxs ## Read the data as a tibble using the format specified at ## https://github.com/ChristopherWilks/megadepth#megadepth-pathtobamfile---junctions example_jxs <- read_junction_table(example_jxs) example_jxs process_junction_table(example_jxs) ## ----createVignette, eval=FALSE----------------------------------------------- # ## Create the vignette # library("rmarkdown") # system.time(render("megadepth.Rmd", "BiocStyle::html_document")) # # ## Extract the R code # library("knitr") # knit("megadepth.Rmd", tangle = TRUE) ## ----reproduce1, echo=FALSE--------------------------------------------------- ## Date the vignette was generated Sys.time() ## ----reproduce2, echo=FALSE--------------------------------------------------- ## Processing time in seconds totalTime <- diff(c(startTime, Sys.time())) round(totalTime, digits = 3) ## ----reproduce3, echo=FALSE------------------------------------------------------------------------------------------- ## Session info library("sessioninfo") options(width = 120) session_info() ## ----vignetteBiblio, results = "asis", echo = FALSE, warning = FALSE, message = FALSE--------------------------------- ## Print bibliography PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))