getStrandedCoverage {BRGenomics} | R Documentation

## Get strand-specific coverage

### Description

Computes strand-specific coverage signal, and returns a GRanges object. Function also works for non-strand-specific data.

### Usage

    getStrandedCoverage(
      dataset.gr,
      field = "score",
      ncores = getOption("mc.cores", 2L)
    )

### Arguments

`dataset.gr` | A GRanges object either containing ranges for each read, or one in which readcounts for individual ranges are contained in metadata (typically in the "score" field). `dataset.gr` can also be a list of such GRanges objects.
`field` | The name of the metadata field that contains readcounts. If no metadata field contains readcounts, and each range represents a single read, set to NULL.
`ncores` | Number of cores to use for calculating coverage. For a single dataset, the max that will be used is 3, one for each possible strand (plus, minus, and unstranded). More cores can be used if `dataset.gr` is a list.

### Value

A GRanges object with signal in the "score" metadata column. Note that the output is not automatically converted into a `"basepair-resolution"` GRanges object.

### Author(s)

Mike DeBerardine

### See Also

`makeGRangesBRG`, `GenomicRanges::coverage`

### Examples

    #--------------------------------------------------#
    # Using included full-read data
    #--------------------------------------------------#
    # -> whole-read coverage sacrifices meaningful readcount
    #    information, but can be useful for visualization,
    #    e.g. for looking at RNA-seq data in a genome browser

    data("PROseq_paired")

    PROseq_paired[1:6]

    getStrandedCoverage(PROseq_paired, ncores = 1)[1:6]

    #--------------------------------------------------#
    # Getting coverage from single bases of single reads
    #--------------------------------------------------#

    # included PROseq data is already single-base coverage
    data("PROseq")
    range(width(PROseq))

    # undo coverage for the first 100 positions
    ps <- PROseq[1:100]
    ps_reads <- rep(ps, times = ps$score)
    mcols(ps_reads) <- NULL

    ps_reads[1:6]

    # re-create coverage
    getStrandedCoverage(ps_reads, field = NULL, ncores = 1)[1:6]

    #--------------------------------------------------#
    # Reversing makeGRangesBRG
    #--------------------------------------------------#
    # -> getStrandedCoverage doesn't return single-width
    #    GRanges, which is useful because getting coverage
    #    will merge adjacent bases with equivalent scores

    # included PROseq data is already single-width
    range(width(PROseq))
    isDisjoint(PROseq)

    ps_cov <- getStrandedCoverage(PROseq, ncores = 1)

    range(width(ps_cov))
    sum(score(PROseq)) == sum(score(ps_cov) * width(ps_cov))

    # -> Look specifically at ranges that could be combined
    neighbors <- c(shift(PROseq, 1), shift(PROseq, -1))
    hits <- findOverlaps(PROseq, neighbors)
    idx <- unique(from(hits)) # indices for PROseq with neighbor

    PROseq[idx]

    getStrandedCoverage(PROseq[idx], ncores = 1)

[Package BRGenomics version 1.2.0 Index]