GappedAlignments-class {GenomicRanges}R Documentation

GappedAlignments objects

Description

The GappedAlignments class is a simple container which purpose is to store a set of alignments that will hold just enough information for supporting the operations described below.

Details

A GappedAlignments object is a vector-like object where each element describes an alignment i.e. how a given sequence (called "query" or "read", typically short) aligns to a reference sequence (typically long).

Typically, a GappedAlignments object will be created by loading records from a BAM (or SAM) file and each element in the resulting object will correspond to a record. BAM/SAM records generally contain a lot of information but only part of that information is loaded in the GappedAlignments object. In particular, we discard the query sequences (SEQ field), the query qualities (QUAL), the mapping qualities (MAPQ) and any other information that is not needed in order to support the operations or methods described below.

This means that multi-reads (i.e. reads with multiple hits in the reference) won't receive any special treatment i.e. the various SAM/BAM records corresponding to a multi-read will show up in the GappedAlignments object as if they were coming from different/unrelated queries. Also paired-end reads will be treated as single-end reads and the pairing information will be lost (see ?GappedAlignmentPairs for how to handle aligned paired-end reads).

Each element of a GappedAlignments object consists of:

Note that the last 2 items are not expicitly stored in the GappedAlignments object: they are inferred on-the-fly from the CIGAR and the "start".

Optionally, a GappedAlignments object can have names (accessed thru the names generic function) which will be coming from the QNAME field of the SAM/BAM records.

The rest of this man page will focus on describing how to:

Constructors

readGappedAlignments(file, format="BAM", use.names=FALSE, ...): Read a file containing aligned reads as a GappedAlignments object. By default (i.e. use.names=FALSE), the resulting object has no names. If use.names is TRUE, then the names are constructed from the query template names (QNAME field in a SAM/BAM file).

Note that this function is just a front-end that delegates to the format-specific back-end function specified via the format argument. The use.names argument and any extra argument are passed to the back-end function. Only the BAM format is supported for now. Its back-end is the readBamGappedAlignments function defined in the Rsamtools package. See ?readBamGappedAlignments for more information (you might need to install and load the Rsamtools package first).

GappedAlignments(seqnames=Rle(factor()), pos=integer(0), cigar=character(0), strand=NULL, names=NULL, seqlengths=NULL, ...): Low-level GappedAlignments constructor. Generally not used directly. Named arguments in ... are used as metadata columns.

Accessors

In the code snippets below, x is a GappedAlignments object.

length(x): Return the number of alignments in x.

names(x), names(x) <- value: Get or set the names of x. See readGappedAlignments above for how to automatically extract and set the names from the file to read.

seqnames(x), seqnames(x) <- value: Get or set the name of the reference sequence for each alignment in x (see Details section above for more information about the RNAME field of a SAM/BAM file). value can be a factor, or a 'factor' Rle, or a character vector.

rname(x), rname(x) <- value: Same as seqnames(x) and seqnames(x) <- value.

strand(x), strand(x) <- value: Get or set the strand for each alignment in x (see Details section above for more information about the strand of an alignment). value can be a factor (with levels +, - and *), or a 'factor' Rle, or a character vector.

cigar(x): Returns a character vector of length length(x) containing the CIGAR string for each alignment.

qwidth(x): Returns an integer vector of length length(x) containing the length of the query *after* hard clipping (i.e. the length of the query sequence that is stored in the corresponding SAM/BAM record).

start(x), end(x): Returns an integer vector of length length(x) containing the "start" and "end" (respectively) of the query for each alignment. See Details section above for the exact definitions of the "start" and "end" of a query. Note that start(x) and end(x) are equivalent to start(granges(x)) and end(granges(x)), respectively (or, alternatively, to min(rglist(x)) and max(rglist(x)), respectively).

width(x): Equivalent to width(granges(x)) (or, alternatively, to end(x) - start(x) + 1L). Note that this is generally different from qwidth(x) except for alignments with a trivial CIGAR string (i.e. a string of the form "<n>M" where <n> is a number).

ngap(x): Returns an integer vector of the same length as x containing the number of gaps (i.e. N operations in the CIGAR) for each alignment. Equivalent to unname(elementLengths(rglist(x))) - 1L.

seqinfo(x), seqinfo(x) <- value: Get or set the information about the underlying sequences. value must be a Seqinfo object.

seqlevels(x), seqlevels(x) <- value: Get or set the sequence levels. seqlevels(x) is equivalent to seqlevels(seqinfo(x)) or to levels(seqnames(x)), those 2 expressions being guaranteed to return identical character vectors on a GappedAlignments object. value must be a character vector with no NAs. See ?seqlevels for more information.

seqlengths(x), seqlengths(x) <- value: Get or set the sequence lengths. seqlengths(x) is equivalent to seqlengths(seqinfo(x)). value can be a named non-negative integer or numeric vector eventually with NAs.

isCircular(x), isCircular(x) <- value: Get or set the circularity flags. isCircular(x) is equivalent to isCircular(seqinfo(x)). value must be a named logical vector eventually with NAs.

genome(x), genome(x) <- value: Get or set the genome identifier or assembly name for each sequence. genome(x) is equivalent to genome(seqinfo(x)). value must be a named character vector eventually with NAs.

seqnameStyle(x): Get or set the seqname style for x. Note that this information is not stored in x but inferred by looking up seqnames(x) against a seqname style database stored in the seqnames.db metadata package (required). seqnameStyle(x) is equivalent to seqnameStyle(seqinfo(x)) and can return more than 1 seqname style (with a warning) in case the style cannot be determined unambiguously.

Coercion

In the code snippets below, x is a GappedAlignments object.

grglist(x, order.as.in.query=FALSE, drop.D.ranges=FALSE), rglist(x, order.as.in.query=FALSE, drop.D.ranges=FALSE):

Return either a GRangesList or a RangesList object of length length(x) where the i-th element represents the ranges (with respect to the reference) of the i-th alignment in x.

More precisely, the RangesList object returned by rglist(x) is a CompressedIRangesList object.

The order.as.in.query toggle affects the order of the ranges within each top-level element of the returned object.

If FALSE (the default), then the ranges are ordered from 5' to 3' in elements associated with the plus strand (i.e. corresponding to alignments located on the plus strand), and from 3' to 5' in elements associated with the minus strand. So, whatever the strand is, the ranges are in ascending order (i.e. left-to-right).

If TRUE, then the order of the ranges in elements associated with the minus strand is reversed. So they end up being ordered from 5' to 3' too, which means that they are now in decending order (i.e. right-to-left). It also means that, when order.as.in.query=TRUE is used, the ranges are always ordered consistently with the original "query template", that is, in the order defined by walking the "query template" from the beginning to the end.

If drop.D.ranges is TRUE, then deletions (D operations in the CIGAR) are treated like gaps (N operations in the CIGAR), that is, the ranges corresponding to deletions are dropped.

See Details section above for more information.

granges(x), ranges(x): Return either a GRanges or a Ranges object of length length(x) where each element represents the regions in the reference to which a query is aligned.

More precisely, the Ranges object returned by ranges(x) is an IRanges object.

introns(x): Extract the gaps (i.e. N operations in the CIGAR) as a GRangesList object of the same length as x. Equivalent to:

    psetdiff(granges(x), grglist(x, order.as.in.query=TRUE))
      

as(x, "GRangesList"), as(x, "GRanges"), as(x, "RangesList"), as(x, "Ranges"): Alternate ways of doing grglist(x), granges(x), rglist(x), ranges(x), respectively.

Subsetting and related operations

In the code snippets below, x is a GappedAlignments object.

x[i]: Return a new GappedAlignments object made of the selected alignments. i can be a numeric or logical vector.

Combining

c(...): Concatenates the GappedAlignment objects in ....

Other methods

qnarrow(x, start=NA, end=NA, width=NA): x is a GappedAlignments object. Return a new GappedAlignments object of the same length as x describing how the narrowed query sequences align to the reference. The start/end/width arguments describe how to narrow the query sequences. They must be vectors of integers. NAs and negative values are accepted and "solved" according to the rules of the SEW (Start/End/Width) interface (see ?solveUserSEW for the details).

show(x): By default the show method displays 5 head and 5 tail lines. The number of lines can be altered by setting the global options showHeadLines and showTailLines. If the object length is less than the sum of the options, the full object is displayed. These options affect GRanges, GappedAlignments, Ranges, DataTable and XString objects.

Author(s)

H. Pages and P. Aboyoun

References

http://samtools.sourceforge.net/

See Also

Examples

library(Rsamtools)  # for ScanBamParam() and the ex1.bam file
ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
gal <- readGappedAlignments(ex1_file, param=ScanBamParam(what="flag"))
gal

## ---------------------------------------------------------------------
## A. BASIC MANIPULATION
## ---------------------------------------------------------------------
length(gal)
head(gal)
names(gal)  # no names by default
seqnames(gal)
strand(gal)
head(cigar(gal))
head(qwidth(gal))
table(qwidth(gal))
head(start(gal))
head(end(gal))
head(width(gal))
head(ngap(gal))
seqlevels(gal)

## Rename the reference sequences:
seqlevels(gal) <- sub("seq", "chr", seqlevels(gal))
seqlevels(gal)

grglist(gal)  # a GRangesList object
stopifnot(identical(unname(elementLengths(grglist(gal))), ngap(gal) + 1L))
granges(gal)  # a GRanges object
rglist(gal)   # a CompressedIRangesList object
stopifnot(identical(unname(elementLengths(rglist(gal))), ngap(gal) + 1L))
ranges(gal)   # an IRanges object
introns(gal)  # a GRangesList object
stopifnot(identical(unname(elementLengths(introns(gal))), ngap(gal)))

## Modify the number of lines in 'show'
options("showHeadLines"=3)
options("showTailLines"=2)
gal

## Revert to default
options("showHeadLines"=NULL)
options("showTailLines"=NULL)

## ---------------------------------------------------------------------
## B. SUBSETTING
## ---------------------------------------------------------------------
gal[strand(gal) == "-"]
gal[grep("I", cigar(gal), fixed=TRUE)]
gal[grep("N", cigar(gal), fixed=TRUE)]  # no gaps

## A confirmation that all the queries map to the reference with no
## gaps:
stopifnot(all(ngap(gal) == 0))

## Different ways to subset:
gal[6]             # a GappedAlignments object of length 1
grglist(gal)[[6]]  # a GRanges object of length 1
rglist(gal)[[6]]   # a NormalIRanges object of length 1

## D operations are NOT gaps:
ii <- grep("D", cigar(gal), fixed=TRUE)
gal[ii]
ngap(gal[ii])
grglist(gal[ii])

## qwidth() vs width():
gal[qwidth(gal) != width(gal)]

## This MUST return an empty object:
gal[cigar(gal) == "35M" & qwidth(gal) != 35]
## but this doesn't have too:
gal[cigar(gal) != "35M" & qwidth(gal) == 35]

## ---------------------------------------------------------------------
## C. qnarrow()/narrow()
## ---------------------------------------------------------------------
## Note that there is no difference between qnarrow() and narrow() when
## all the alignments are simple and with no indels.

## This trims 3 nucleotides on the left and 5 nucleotides on the right
## of each alignment:
qnarrow(gal, start=4, end=-6)
## Note that the 'start' and 'end' arguments specify what part of each
## query sequence should be kept (negative values being relative to the
## right end of the query sequence), not what part should be trimmed.

## Trimming on the left doesn't change the "end" of the queries.
qnarrow(gal, start=21)
stopifnot(identical(end(qnarrow(gal, start=21)), end(gal)))

[Package GenomicRanges version 1.12.1 Index]