get_coverage: Compute coverage summarizations across a set of regions
In LieberInstitute/megadepth: megadepth: BigWig and BAM related utilities

View source: R/get_coverage.R

get_coverage

R Documentation

Compute coverage summarizations across a set of regions

Description

Given an input set of annotation regions, compute coverage summarizations using Megadepth for a given BigWig file.

Usage

get_coverage(
  bigwig_file,
  op = c("sum", "mean", "max", "min"),
  annotation,
  prefix = file.path(tempdir(), "bw.mean")
)

Arguments

`bigwig_file`	A `character(1)` with the path to the input BigWig file.
`op`	A `character(1)` specifying the summarization operation to perform.
`annotation`	A `character(1)` path to a BED file with the genomic coordinates you are interested in.
`prefix`	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.

Details

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.

Value

A GRanges-class object with the coverage summarization across the annotation ranges.

Examples


## 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,
    )
}

LieberInstitute/megadepth documentation built on May 6, 2024, 7:12 p.m.