Nothing
## ----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 <https://github.com/LieberInstitute/megadepth/tree/master/analysis> 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
read_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"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.