coverage-methods {GenomicAlignments} | R Documentation

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

### Usage

    ## S4 method for signature 'GAlignments'
    coverage(x, shift=0L, width=NULL, weight=1L,
             method=c("auto", "sort", "hash"), drop.D.ranges=FALSE)

    ## S4 method for signature 'GAlignmentPairs'
    coverage(x, shift=0L, width=NULL, weight=1L,
             method=c("auto", "sort", "hash"), 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)

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

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

### Value

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

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

[Package GenomicAlignments version 1.20.1 Index]