R/read_coverage.R

In megadepth: megadepth: BigWig and BAM related utilities

#' Read a coverage TSV file created by Megadepth
#'
#' Read an `*annotation.tsv` file created by `get_coverage()` or manually by
#' the user using Megadepth.
#'
#' @param tsv_file A `character(1)` specifying the path to the tab-separated
#' (TSV) file created manually using `megadepth_shell()` or on a previous
#' `get_coverage()` run.
#' @param verbose A `logical(1)` controlling whether to suppress messages when
#' reading the data.
#'
#' @return A [GRanges-class][GenomicRanges::GRanges-class] object with the
#' coverage summarization across the annotation ranges.
#' @export
#' @family Coverage functions
#' @importFrom GenomicRanges GRanges
#' @importFrom readr read_delim cols
#'
#' @examples
#'
#' ## Install if necessary
#' install_megadepth()
#'
#' ## Locate 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
#'
#' ## Read in the coverage file again, using read_coverage()
#' ## First, lets locate the tsv file that was generated by get_coverage()
#' tsv_file <- file.path(tempdir(), "bw.mean.annotation.tsv")
#' bw_cov_manual <- read_coverage(tsv_file)
#' stopifnot(identical(bw_cov, bw_cov_manual))
#'
#' ## To get an RleList object, just like the one you would get
#' ## from using rtracklayer::import.bw(as = "RleList") directly on the
#' ## BigWig file, use:
#' GenomicRanges::coverage(bw_cov_manual)
#'
#' ## The coverage data can also be read as a `tibble::tibble()`
#' read_coverage_table(tsv_file)
read_coverage <- function(tsv_file, verbose = TRUE) {
    if (verbose) {
        coverage <- read_coverage_table(tsv_file)
    } else {
        suppressMessages(coverage <- read_coverage_table(tsv_file))
    }

    ## Cast into a GRanges object
    coverage <- GenomicRanges::GRanges(coverage)

    return(coverage)
}


#' @describeIn read_coverage Read a coverage TSV file created by Megadepth as
#' a table
#'
#' @return A `tibble::tible()` with columns `chr`, `start`, `end` and `score`.
#' @export
read_coverage_table <- function(tsv_file) {
    coverage <- readr::read_delim(
        tsv_file,
        delim = "\t",
        col_names = c("chr", "start", "end", "score"),
        col_types = readr::cols(
            chr = "c",
            start = "i",
            end = "i",
            score = "n"
        ),
        progress = FALSE
    )
    return(coverage)
}