get_coverage {megadepth} | R Documentation |
Given an input set of annotation regions, compute coverage summarizations using Megadepth for a given BigWig file.
get_coverage( bigwig_file, op = c("sum", "mean", "max", "min"), annotation, prefix = file.path(tempdir(), "bw.mean") )
bigwig_file |
A |
op |
A |
annotation |
A |
prefix |
A |
Note that the chromosome names (seqnames) in the BigWig file and the annotation file should use the same format. Otherwise, Megadepth will return 0 counts.
A GRanges-class object with the coverage summarization across the annotation ranges.
Other Coverage functions:
read_coverage()
## Install if necessary install_megadepth() ## 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 ) ## Compute the coverage bw_cov <- get_coverage(example_bw, op = "mean", annotation = annotation_file) bw_cov ## If you want to cast this into a RleList object use the following code: ## (it's equivalent to rtracklayer::import.bw(as = "RleList")) ## although in the megadepth case the data has been summarized GenomicRanges::coverage(bw_cov) ## 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 if (!xfun::is_windows()) { 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, ) }