Coverage of a GAlignments, GAlignmentPairs, or GAlignmentsList object
coverage
methods for GAlignments,
GAlignmentPairs, GAlignmentsList, and
BamFile objects.
NOTE: The coverage
generic function and methods
for IntegerRanges and IntegerRangesList
objects are defined and documented in the IRanges package.
Methods for GRanges and
GRangesList objects are defined and
documented in the GenomicRanges package.
## S4 method for signature 'GAlignments' coverage(x, shift=0L, width=NULL, weight=1L, method=c("auto", "sort", "hash", "naive"), drop.D.ranges=FALSE) ## S4 method for signature 'GAlignmentPairs' coverage(x, shift=0L, width=NULL, weight=1L, method=c("auto", "sort", "hash", "naive"), drop.D.ranges=FALSE) ## S4 method for signature 'GAlignmentsList' coverage(x, shift=0L, width=NULL, weight=1L, ...) ## S4 method for signature 'BamFile' coverage(x, shift=0L, width=NULL, weight=1L, ..., param=ScanBamParam()) ## S4 method for signature 'character' coverage(x, shift=0L, width=NULL, weight=1L, ..., yieldSize=2500000L)
x |
A GAlignments, GAlignmentPairs, GAlignmentsList, or BamFile object, or the path to a BAM file. |
shift, width, weight |
See |
method |
See |
drop.D.ranges |
Whether the coverage calculation should ignore ranges corresponding to D (deletion) in the CIGAR string. |
... |
Additional arguments passed to the |
param |
An optional ScanBamParam object passed to
|
yieldSize |
An optional argument controlling how many records are input when iterating through a BamFile. |
The methods for GAlignments and GAlignmentPairs objects do:
coverage(grglist(x, drop.D.ranges=drop.D.ranges), ...)
The method for GAlignmentsList objects does:
coverage(unlist(x), ...)
The method for BamFile objects iterates through a
BAM file, reading yieldSize(x)
records (or all records, if
is.na(yieldSize(x))
) and calculating:
gal <- readGAlignments(x, param=param) coverage(gal, shift=shift, width=width, weight=weight, ...)
The method for character
vectors of length 1 creates a
BamFile object from x
and performs the
calculation for coverage,BamFile-method
.
A named RleList object with one coverage vector per
seqlevel in x
.
coverage
in the IRanges package.
coverage-methods in the GenomicRanges package.
RleList objects in the IRanges package.
GAlignments and GAlignmentPairs objects.
BamFile objects in the Rsamtools package.
## --------------------------------------------------------------------- ## A. EXAMPLE WITH TOY DATA ## --------------------------------------------------------------------- ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools") ## Coverage of a GAlignments object: gal <- readGAlignments(ex1_file) cvg1 <- coverage(gal) cvg1 ## Coverage of a GAlignmentPairs object: galp <- readGAlignmentPairs(ex1_file) cvg2 <- coverage(galp) cvg2 ## Coverage of a GAlignmentsList object: galist <- readGAlignmentsList(ex1_file) cvg3 <- coverage(galist) cvg3 table(mcols(galist)$mate_status) mated_idx <- which(mcols(galist)$mate_status == "mated") mated_galist <- galist[mated_idx] mated_cvg3 <- coverage(mated_galist) mated_cvg3 ## Sanity checks: stopifnot(identical(cvg1, cvg3)) stopifnot(identical( cvg2, mated_cvg3)) ## --------------------------------------------------------------------- ## B. EXAMPLE WITH REAL DATA ## --------------------------------------------------------------------- library(pasillaBamSubset) ## See '?pasillaBamSubset' for more information about the 2 BAM files ## included in this package. reads <- readGAlignments(untreated3_chr4()) table(njunc(reads)) # data contains junction reads ## Junctions do NOT contribute to the coverage: read1 <- reads[which(njunc(reads) != 0L)[1]] # 1st read with a junction read1 # cigar shows a "skipped region" of length 15306 grglist(read1)[[1]] # the junction is between pos 4500 and 19807 coverage(read1)$chr4 # junction is not covered ## Sanity checks: cvg <- coverage(reads) read_chunks <- unlist(grglist(reads), use.names=FALSE) read_chunks_per_chrom <- split(read_chunks, seqnames(read_chunks)) stopifnot(identical(sum(cvg), sum(width(read_chunks_per_chrom)))) galist <- readGAlignmentsList(untreated3_chr4()) stopifnot(identical(cvg, coverage(galist)))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.