inst/doc/megadepth.R

In megadepth: megadepth: BigWig and BAM related utilities

## ----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"))