summarizeOverlaps {GenomicRanges}R Documentation

Perform overlap queries between reads and genomic features

Description

summarizeOverlaps extends findOverlaps by providing options to resolve reads that overlap multiple features.

Usage

  ## S4 method for signature 'GRanges,GAlignments'
summarizeOverlaps(
    features, reads, mode, ignore.strand=FALSE, ..., inter.feature=TRUE) 
  ## S4 method for signature 'GRangesList,GAlignments'
summarizeOverlaps(
    features, reads, mode, ignore.strand=FALSE, ..., inter.feature=TRUE) 
  ## S4 method for signature 'GRanges,GAlignmentPairs'
summarizeOverlaps(
    features, reads, mode, ignore.strand=FALSE, ..., inter.feature=TRUE) 
  ## S4 method for signature 'GRangesList,GAlignmentPairs'
summarizeOverlaps(
    features, reads, mode, ignore.strand=FALSE, ..., inter.feature=TRUE) 

  ## mode funtions
  Union(features, reads, ignore.strand=FALSE, inter.feature=TRUE)
  IntersectionStrict(features, reads, ignore.strand=FALSE, inter.feature=TRUE)
  IntersectionNotEmpty(features, reads, ignore.strand=FALSE, inter.feature=TRUE)

Arguments

reads

A BamFileList that represents the data to be counted by summarizeOverlaps.

features

A GRanges or a GRangesList object of genomic regions of interest. When a GRanges is supplied, each row is considered a feature. When a GRangesList is supplied, each higher list-level is considered a feature. This distinction is important when defining overlaps.

mode

A function that defines the method to be used when a read overlaps more than one feature. Pre-defined options are "Union", "IntersectionStrict", or "IntersectionNotEmpty" and are designed after the counting modes available in the HTSeq package by Simon Anders (see references).

  • "Union" : (Default) Reads that overlap any portion of exactly one feature are counted. Reads that overlap multiple features are discarded. This is the most conservative of the 3 modes.

  • "IntersectionStrict" : A read must fall completely "within" the feature to be counted. If a read overlaps multiple features but falls "within" only one, the read is counted for that feature. If the read is "within" multiple features, the read is discarded.

  • "IntersectionNotEmpty" : A read must fall in a unique disjoint region of a feature to be counted. When a read overlaps multiple features, the features are partitioned into disjoint intervals. Regions that are shared between the features are discarded leaving only the unique disjoint regions. If the read overlaps one of these remaining regions, it is assigned to the feature the unique disjoint region came from.

  • user supplied function : A function can be supplied as the mode argument. It must (1) have arguments that correspond to features, reads, ignore.strand and inter.feature arguments (as in the defined mode functions) and (2) return a vector of counts the same length as features.

ignore.strand

A logical indicating if strand should be considered when matching.

inter.feature

(Default TRUE) A logical indicating if the counting mode should be aware of overlapping features. When TRUE (default), reads mapping to multiple features are dropped (i.e., not counted). When FALSE, these reads are retained and a count is assigned to each feature they map to.

There are 6 possible combinations of the mode and inter.feature arguments. When inter.feature=FALSE the behavior of modes ‘Union’ and ‘IntersectionStrict’ are essentially ‘countOverlaps’ with ‘type=any’ and type=within, respectively. ‘IntersectionNotEmpty’ does not reduce to a simple countOverlaps because common (shared) regions of the annotation are removed before counting.

...

Additional arguments for Bam file methods such as fragments, singleEnd or param. If using multiple cores, arguments can be passed through to mclapply used when counting Bam files.

fragments

(Default FALSE) Applies to paired-end data only so singleEnd must be FALSE.

fragments is a logical value indicating if singletons, reads with unmapped pairs and other fragments should be included in counting. When fragments=FALSE readGAlignmentPairs is used to read in the data, when fragments=TRUE readGAlignmentsList is used.

readGAlignmentPairs keeps only the read pairs mated by the algorithm while readGAlignmentsList keeps the pairs as well as all singletons, reads with unmapped pairs and other fragments. When fragments=TRUE counts will generally be higher because all records are included in the counting, not just the primary alignment pairs. See ?readGAlignmentsListFromBam for the algorithm details.

singleEnd

(Default TRUE) A logical value indicating if reads are single or paired-end. In Bioconductor > 2.12 it is not necessary to sort paired-end Bam files by qname. When counting with summarizeOverlaps, setting singleEnd=FALSE will trigger paired-end reading and counting. It is fine to also set asMates=TRUE in the BamFile but is not necessary when singleEnd=FALSE.

param

An optional ScanBamParam instance to further influence scanning, counting, or filtering.

Details

summarizeOverlaps offers counting modes to resolve reads that overlap multiple features. The mode argument defines a set of rules to resolve the read to a single feature such that each read is counted a maximum of once. New to GenomicRanges >= 1.13.9 is the inter.feature argument which allows reads to be counted for each feature they overlap. When inter.feature=TRUE the counting modes are aware of feature overlap and reads overlapping multiple features are dropped and not counted. When inter.feature=FALSE multiple feature overlap is ignored and reads are counted once for each feature they map to. This essentially reduces modes ‘Union’ and ‘IntersectionStrict’ to countOverlaps with type="any", and type="within", respectively. ‘IntersectionNotEmpty’ is not reduced to a derivative of countOverlaps because the shared regions are removed before counting.

features :

A ‘feature’ can be any portion of a genomic region such as a gene, transcript, exon etc. When the features argument is a GRanges the rows define the features. The result will be the same length as the GRanges. When features is a GRangesList the highest list-level defines the features and the result will be the same length as the GRangesList.

When inter.feature=TRUE, each count mode attempts to assign a read that overlaps multiple features to a single feature. If there are ranges that should be considered together (e.g., exons by transcript or cds regions by gene) the GRangesList would be appropriate. If there is no grouping in the data then a GRanges would be appropriate.

paired-end reads :

Paired-end reads are counted the same as single-end reads with gaps; each pair registers as a single hit. Paired-end records can be counted in a GAlignmentPairs container or Bam file.

When counting Bam files the method has an additional argument, singleEnd, which should be FALSE for paired-end data. See ?summarizeOverlaps,GRanges,BamFileList-method for details on additional arguments when counting Bam files.

Value

A SummarizedExperiment object. The assays slot holds the counts, rowData holds the annotation specified in features.

colData is a DataFrame with columns of ‘object’ (class of reads) and ‘records’ (length of reads). When reads is a BamFile or BamFileList the colData holds the output of a call to countBam with columns of ‘records’ (total records in file), ‘nucleotides’ and ‘mapped’. The number in ‘mapped’ is the number of records returned when isUnmappedQuery=FALSE in the ‘ScanBamParam’.

Author(s)

Valerie Obenchain <vobencha@fhcrc.org>

References

HTSeq : http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html

htseq-count : http://www-huber.embl.de/users/anders/HTSeq/doc/count.html

See Also

Examples

reads <- GAlignments(
    names = c("a","b","c","d","e","f","g"),
    seqnames = Rle(c(rep(c("chr1", "chr2"), 3), "chr1")),
    pos = as.integer(c(1400, 2700, 3400, 7100, 4000, 3100, 5200)),
    cigar = c("500M", "100M", "300M", "500M", "300M", 
              "50M200N50M", "50M150N50M"),
    strand = strand(rep("+", 7)))

gr <- GRanges(
    seqnames = c(rep("chr1", 7), rep("chr2", 4)), strand = "+", 
    ranges = IRanges(c(1000, 3000, 3600, 4000, 4000, 5000, 5400, 
                       2000, 3000, 7000, 7500), 
                     width = c(500, 500, 300, 500, 900, 500, 500, 
                               900, 500, 600, 300),
                     names=c("A", "B", "C1", "C2", "D1", "D2", "E", "F",
                             "G", "H1", "H2"))) 
groups <- factor(c(1,2,3,3,4,4,5,6,7,8,8))
grl <- splitAsList(gr, groups)
names(grl) <- LETTERS[seq_along(grl)]

## ---------------------------------------------------------------------
## Counting modes. 
## ---------------------------------------------------------------------

## First we count with a GRanges as the 'features'. Note that
## 'Union' is the most conservative counting mode followed by 
## 'IntersectionStrict' then 'IntersectionNotEmpty'.
counts1 <- 
    data.frame(union=assays(summarizeOverlaps(gr, reads))$counts, 
               intStrict=assays(summarizeOverlaps(gr, reads, 
                                mode="IntersectionStrict"))$counts,
               intNotEmpty=assays(summarizeOverlaps(gr, reads,
                                  mode="IntersectionNotEmpty"))$counts)

colSums(counts1)

## Split the 'features' into a GRangesList and count again.
counts2 <- 
    data.frame(union=assays(summarizeOverlaps(grl, reads))$counts, 
               intStrict=assays(summarizeOverlaps(grl, reads, 
                                mode="IntersectionStrict"))$counts,
               intNotEmpty=assays(summarizeOverlaps(grl, reads,
                                  mode="IntersectionNotEmpty"))$counts)
colSums(counts2)

## The GRangesList ('grl' object) has 8 features whereas the GRanges 
## ('gr' object) has 11. The affect on counting can be seen by looking
## at feature 'H' with mode 'Union'. In the GRanges this feature is 
## represented by ranges 'H1' and 'H2',
gr[c("H1", "H2")]

## and by list element 'H' in the GRangesList, 
grl["H"]
 
## Read "d" hits both 'H1' and 'H2'. This is considered a multi-hit when
## using a GRanges (each range is a separate feature) so the read was 
## dropped and not counted.
counts1[c("H1", "H2"), ]

## When using a GRangesList, each list element is considered a feature.
## The read hits multiple ranges within list element 'H' but only one 
## list element. This is not considered a multi-hit so the read is counted.
counts2["H", ]

## ---------------------------------------------------------------------
## Counting multi-hit reads.
## ---------------------------------------------------------------------

## The goal of the counting modes is to provide a set of rules that
## resolve reads hitting multiple features so each read is counted
## a maximum of once. However, sometimes it may be desirable to count 
## a read for each feature it overlaps. This can be accomplished by 
## setting 'inter.feature' to FALSE.

## When 'inter.feature=FALSE', modes 'Union' and 'IntersectionStrict'
## essentially reduce to countOverlaps() with type="any" and 
## type="within", respectively.

## When 'inter.feature=TRUE' only features "A", "F" and "G" have counts.
se1 <- summarizeOverlaps(gr, reads, mode="Union", inter.feature=TRUE)
assays(se1)$counts

## When 'inter.feature=FALSE' all 11 features have a count. There are 
## 7 total reads so one or more reads were counted more than once.
se2 <- summarizeOverlaps(gr, reads, mode="Union", inter.feature=FALSE)
assays(se2)$counts

## ---------------------------------------------------------------------
## Counting Bam files.
## ---------------------------------------------------------------------

library(Rsamtools)
library(pasillaBamSubset)
library("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
exbygene <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, "gene")

## (i) Single-end :

## Large files can be iterated over in chunks by setting a
## 'yieldSize' on the BamFile.
bf_s <- BamFile(untreated1_chr4(), yieldSize=50000)
se_s <- summarizeOverlaps(exbygene, bf_s, singleEnd=TRUE)
table(assays(se_s)$counts > 0)

## When a character (file name) is provided as 'reads' instead 
## of a BamFile object summarizeOverlaps() will create a BamFile
## and set a reasonable default 'yieldSize'.

## (ii) Paired-end :

## A paired-end file may contain singletons, reads with unmapped
## pairs or reads with more than two fragments. When 'fragments=FALSE'
## only reads paired by the algorithm are included in the counting. 
nofrag <- summarizeOverlaps(exbygene, untreated3_chr4(), 
                            singleEnd=FALSE, fragments=FALSE)
table(assays(nofrag)$counts > 0)

## When 'fragments=TRUE' all singletons, reads with unmapped pairs 
## and other fragments will be included in the counting.
bf <- BamFile(untreated3_chr4())
frag <- summarizeOverlaps(exbygene, bf, singleEnd=FALSE, fragments=TRUE)
table(assays(frag)$counts > 0)

## As expected, using 'fragments=TRUE' results in a larger number 
## of total counts because singletons, unmapped pairs etc. are 
## included in the counting.

## Total reads in the file:
countBam(untreated3_chr4())

## Reads counted with 'fragments=FALSE':
sum(assays(nofrag)$counts)

## Reads counted with 'fragments=TRUE':
sum(assays(frag)$counts)

## ---------------------------------------------------------------------
## Count tables for DESeq or edgeR.
## ---------------------------------------------------------------------

fls <- list.files(system.file("extdata",package="GenomicRanges"),
                  recursive=TRUE, pattern="*bam$", full=TRUE)
names(fls) <- basename(fls)
bf <- BamFileList(fls, index=character(), yieldSize=1000)
genes <- GRanges(
    seqnames = c(rep("chr2L", 4), rep("chr2R", 5), rep("chr3L", 2)),
    ranges = IRanges(c(1000, 3000, 4000, 7000, 2000, 3000, 3600, 
                       4000, 7500, 5000, 5400), 
                     width=c(rep(500, 3), 600, 900, 500, 300, 900, 
                             300, 500, 500))) 
se <- summarizeOverlaps(genes, bf)

## When the reads are Bam files, the 'colData' contains summary 
## information from a call to countBam().
colData(se)

## Create count tables.
library(DESeq)
deseq <- newCountDataSet(assays(se)$counts, rownames(colData(se)))
library(edgeR)
edger <- DGEList(assays(se)$counts, group=rownames(colData(se)))

## ---------------------------------------------------------------------
## User supplied 'mode'. 
## ---------------------------------------------------------------------

## A user defined count function must have the same arguments as 
## the current counting modes.
## Not run: 
counter <- function(x, y,  ignore.strand, inter.feature) {
   ## count ...
}

se <- summarizeOverlaps(gr, reads, mode=counter) 

## End(Not run)
[Package GenomicRanges version 1.14.4 Index]