csaw (version 1.6.1)

regionCounts: Count reads overlapping each region

Description

Count the number of extended reads overlapping pre-specified regions

Usage

regionCounts(bam.files, regions, ext=100, param=readParam())

Arguments

bam.files

a character vector containing paths to sorted and indexed BAM files

regions

a GRanges object containing the regions over which reads are to be counted

ext

an integer scalar or vector, describing the average length of the sequenced fragment in each library

param

a readParam object containing read extraction parameters, or a list of such objects (one for each BAM file)

Value

A RangedSummarizedExperiment object is returned containing one integer matrix. Each entry of the matrix contains the count for each library (column) at each region (row). The coordinates of each region are stored as the rowRanges. The total number of reads, read extension length and param used in each library are in the colData.

Details

This function simply provides a wrapper around countOverlaps for read counting into specified regions. It is provided so as to allow for counting with awareness of the other parameters, e.g., ext, pe. This allows users to coordinate region-based counts with those from windowCounts. Checking that the output totals are the same between the two calls is strongly recommended.

Note that the strandedness of regions will not be considered when computing overlaps. The strandedness of the output rowRanges will depend on the strand(s) from which reads were counted. This is determined by the forward slot in the param object.

See windowCounts for more details on read extension.

See Also

countOverlaps, windowCounts, readParam

Examples

# A low filter is only used here as the examples have very few reads.
bamFiles <- system.file("exdata", c("rep1.bam", "rep2.bam"), package="csaw")
incoming <- GRanges(c('chrA', 'chrA', 'chrB', 'chrC'), 
    IRanges(c(1, 500, 100, 1000), c(200, 1000, 700, 1500)))
regionCounts(bamFiles, regions=incoming)
regionCounts(bamFiles, regions=incoming, param=readParam(restrict="chrB"))

# Loading PE data.
bamFile <- system.file("exdata", "pet.bam", package="csaw")
regionCounts(bamFile, regions=incoming, param=readParam(pe="both"))
regionCounts(bamFile, regions=incoming, param=readParam(max.frag=100, 
	pe="first", restrict="chrA"))