Read an *annotation.tsv file created by get_coverage() or manually by the user using Megadepth.

read_coverage(tsv_file, verbose = TRUE)

read_coverage_table(tsv_file)

Arguments

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.

verbose

A logical(1) controlling whether to suppress messages when reading the data.

Value

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

A tibble::tibble() with columns chr, start, end and score.

Functions

  • read_coverage_table(): Read a coverage TSV file created by Megadepth as a table

See also

Other Coverage functions: get_coverage()

Examples


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

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

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

## The coverage data can also be read as a `tibble::tibble()`
read_coverage_table(tsv_file)
#> # A tibble: 4 × 4
#>   chr          start     end score
#>   <chr>        <int>   <int> <dbl>
#> 1 chr10            0      10  0   
#> 2 chr10      8756697 8756762 15.8 
#> 3 chr10      4359156 4359188  3   
#> 4 GL000219.1  168500  168620  1.26