Given an input set of annotation regions, compute coverage summarizations using Megadepth for a given BigWig file.
A character(1) with the path to the input BigWig file.
A character(1) specifying the summarization operation to
perform.
A character(1) path to a BED file with the genomic
coordinates you are interested in.
A character(1) specifying the output file prefix. This
function creates a file called prefix.annotation.tsv that can be read
again later with read_coverage(). By default the file is created in
the tempdir() and will be deleted after you close your R session.
A GRanges-class object with the coverage summarization across the annotation ranges.
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.
Other Coverage functions:
read_coverage()
## Install if necessary
install_megadepth()
#> The latest megadepth version is 1.2.0
#> This is not an interactive session, therefore megadepth has been installed temporarily to
#> /tmp/RtmpK16vxZ/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
#> GRanges object with 4 ranges and 1 metadata column:
#> seqnames ranges strand | score
#> <Rle> <IRanges> <Rle> | <numeric>
#> [1] chr10 0-10 * | 0.00
#> [2] chr10 8756697-8756762 * | 15.85
#> [3] chr10 4359156-4359188 * | 3.00
#> [4] GL000219.1 168500-168620 * | 1.26
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 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)
#> RleList of length 2
#> $chr10
#> integer-Rle of length 8756762 with 5 runs
#> Lengths: 10 4359145 33 4397508 66
#> Values : 1 0 1 0 1
#>
#> $GL000219.1
#> integer-Rle of length 168620 with 2 runs
#> Lengths: 168499 121
#> Values : 0 1
#>
## 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,
)
}
#> Warning: `expect_equivalent()` was deprecated in the 3rd edition.
#> ℹ Use expect_equal(ignore_attr = TRUE)
#> Warning: `expect_equivalent()` was deprecated in the 3rd edition.
#> ℹ Use expect_equal(ignore_attr = TRUE)