GenomicAlignments (version 1.8.4)

coverage-methods: Coverage of a GAlignments, GAlignmentPairs, or GAlignmentsList object

Description

coverage methods for GAlignments, GAlignmentPairs, GAlignmentsList, and BamFile objects.

NOTE: The coverage generic function and methods for Ranges and RangesList objects are defined and documented in the IRanges package. Methods for GRanges and GRangesList objects are defined and documented in the GenomicRanges package.

Usage

"coverage"(x, shift=0L, width=NULL, weight=1L, method=c("auto", "sort", "hash"), drop.D.ranges=FALSE)
"coverage"(x, shift=0L, width=NULL, weight=1L, method=c("auto", "sort", "hash"), drop.D.ranges=FALSE)
"coverage"(x, shift=0L, width=NULL, weight=1L, ...)
"coverage"(x, shift=0L, width=NULL, weight=1L, ..., param=ScanBamParam())
"coverage"(x, shift=0L, width=NULL, weight=1L, ..., yieldSize=2500000L)

Arguments

x

A GAlignments, GAlignmentPairs, GAlignmentsList, or BamFile object, or the path to a BAM file.

shift, width, weight

See coverage method for GRanges objects in the GenomicRanges package.

method

See ?coverage in the IRanges package for a description of this argument.

drop.D.ranges

Whether the coverage calculation should ignore ranges corresponding to D (deletion) in the CIGAR string.

...

Additional arguments passed to the coverage method for GAlignments objects.

param

An optional ScanBamParam object passed to readGAlignments.

yieldSize

An optional argument controlling how many records are input when iterating through a BamFile.

Value

A named RleList object with one coverage vector per seqlevel in x.

Details

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.

See Also

• coverage in the IRanges package.

• coverage-methods in the GenomicRanges package.

• RleList objects in the IRanges package.

• GAlignments and GAlignmentPairs objects.

• readGAlignments.

• BamFile objects in the Rsamtools package.

Examples

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