getStrandedCoverage {BRGenomics} | R Documentation |
Computes strand-specific coverage signal, and returns a GRanges object. Function also works for non-strand-specific data.
getStrandedCoverage( dataset.gr, field = "score", ncores = getOption("mc.cores", 2L) )
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). |
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 |
A GRanges object with signal in the "score" metadata column. Note
that the output is not automatically converted into a
"basepair-resolution"
GRanges
object.
Mike DeBerardine
makeGRangesBRG
,
GenomicRanges::coverage
#--------------------------------------------------# # 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)