GappedAlignmentPairs-class {GenomicRanges} | R Documentation |
The GappedAlignmentPairs class is a container for "alignment pairs".
A GappedAlignmentPairs object is a list-like object where each element describes an "alignment pair".
An "alignment pair" is made of a "first" and a "last" alignment, and is formally represented by a GappedAlignments object of length 2. It is typically representing a hit of a paired-end read to the reference genome that was used by the aligner. More precisely, in a given pair, the "first" alignment represents the hit of the first end of the read (aka "first segment in the template", using SAM Spec terminology), and the "last" alignment represents the hit of the second end of the read (aka "last segment in the template", using SAM Spec terminology).
In general, a GappedAlignmentPairs object will be created by loading
records from a BAM (or SAM) file containing aligned paired-end reads,
using the readGappedAlignmentPairs
function (see below).
Each element in the returned object will be obtained by pairing 2
records.
readGappedAlignmentPairs(file, format="BAM", use.names=FALSE, ...)
:
Read a file containing paired-end reads as a GappedAlignmentPairs
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 the 2 records in a pair of records have the same QNAME.
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
readBamGappedAlignmentPairs
function
defined in the Rsamtools package.
See ?readBamGappedAlignmentPairs
for
more information (you might need to install and load the Rsamtools
package first).
GappedAlignmentPairs(first, last, isProperPair, names=NULL)
:
Low-level GappedAlignmentPairs constructor. Generally not used directly.
In the code snippets below, x
is a GappedAlignmentPairs object.
length(x)
:
Return the number of alignment pairs in x
.
names(x)
, names(x) <- value
:
Get or set the names of x
.
See readGappedAlignmentPairs
above for how to automatically
extract and set the names from the file to read.
first(x, invert.strand=FALSE)
,
last(x, invert.strand=FALSE)
:
Get the "first" or "last" alignment for each alignment pair in
x
.
The result is a GappedAlignments object of the same length
as x
.
If invert.strand=TRUE
, then the strand is inverted on-the-fly,
i.e. "+" becomes "-", "-" becomes "+", and "*" remains unchanged.
left(x)
:
Get the "left" alignment for each alignment pair in x
.
By definition, the "left" alignment in a pair is the alignment that
is on the + strand. If this is the "first" alignment, then it's returned
as-is by left(x)
, but if this is the "last" alignment, then it's
returned by left(x)
with the strand inverted.
right(x)
:
Get the "right" alignment for each alignment pair in x
.
By definition, the "right" alignment in a pair is the alignment that
is on the - strand. If this is the "first" alignment, then it's returned
as-is by right(x)
, but if this is the "last" alignment, then it's
returned by right(x)
with the strand inverted.
seqnames(x)
:
Get the name of the reference sequence for each alignment pair
in x
. This comes from the RNAME field of the BAM file and
has the same value for the 2 records in a pair
(makeGappedAlignmentPairs
, the function
used by readBamGappedAlignmentPairs
for
doing the pairing, rejects pairs with incompatible RNAME values).
strand(x)
, strand(x) <- value
:
Get or set the strand for each alignment pair in x
.
By definition (and in a somewhat arbitrary way) the strand of an
alignment pair is the strand of the "first" alignment in the pair.
In a GappedAlignmentPairs object, the strand of the "last" alignment
in a pair is typically (but not always) the opposite of the strand
of the "first" alignment. Note that, currently,
makeGappedAlignmentPairs
, the function
used internally by readBamGappedAlignmentPairs
for doing the pairing, rejects pairs where the "first" and "last"
alignments are on the same strand, but those pairs might be supported
in the future.
ngap(x)
:
Equivalent to ngap(first(x)) + ngap(last(x))
.
isProperPair(x)
:
Get the "isProperPair" flag bit (bit 0x2 in SAM Spec) set by
the aligner for each alignment pair in x
.
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
GappedAlignmentPairs 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.
In the code snippets below, x
is a GappedAlignmentPairs object.
x[i]
:
Return a new GappedAlignmentPairs object made of the selected
alignment pairs.
In the code snippets below, x
is a GappedAlignmentPairs object.
x[[i]]
:
Extract the i-th alignment pair as a GappedAlignments object
of length 2. As expected x[[i]][1]
and x[[i]][2]
are
respectively the "first" and "last" alignments in the pair.
unlist(x, use.names=TRUE)
:
Return the GappedAlignments object conceptually defined
by c(x[[1]], x[[2]], ..., x[[length(x)]])
.
use.names
determines whether x
names should be
propagated to the result or not.
In the code snippets below, x
is a GappedAlignmentPairs object.
grglist(x, order.as.in.query=FALSE, drop.D.ranges=FALSE)
:
Return a GRangesList object of length length(x)
where the i-th element represents the ranges (with respect to the
reference) of the i-th alignment pair in x
.
IMPORTANT: The strand of the ranges coming from the "last" alignment in the pair is always inverted.
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 "left" ranges are placed before
the "right" ranges, and, within each left or right group, are ordered
from 5' to 3' in elements associated with the plus strand and from 3'
to 5' in elements associated with the minus strand.
More formally, the i-th element in the returned GRangesList
object can be defined as c(grl1[[i]], grl2[[i]])
, where
grl1
is grglist(left(x))
and grl2
is
grglist(right(x))
.
If TRUE
, then the "first" ranges are placed before the "last"
ranges, and, within each first or last group, are always
ordered from 5' to 3', whatever the strand is.
More formally, the i-th element in the returned GRangesList
object can be defined as c(grl1[[i]], grl2[[i]])
, where
grl1
is grglist(first(x),
order.as.in.query=TRUE)
and
grl2
is grglist(last(x, invert.strand=TRUE),
order.as.in.query=TRUE)
.
Note that the relationship between the 2 GRangesList objects
obtained with order.as.in.query
being respectively
FALSE
or TRUE
is simpler than it sounds: the only
difference is that the order of the ranges in elements associated
with the minus strand is reversed.
Finally note that, in the latter, 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 (Ds in the
CIGAR) are treated like gaps (Ns in the CIGAR), that is, the ranges
corresponding to deletions are dropped.
granges(x)
: Return a GRanges object of length
length(x)
where each range is obtained by merging all the
ranges within the corresponding top-level element in grglist(x)
.
introns(x)
: Extract the gaps (i.e. N operations in the CIGAR)
of the "first" and "last" alignments of each pair as a
GRangesList object of the same length as x
.
Equivalent to (but faster than):
introns1 <- introns(first(x)) introns2 <- introns(last(x, invert.strand=TRUE)) mendoapply(c, introns1, introns2)
as(x, "GRangesList")
, as(x, "GRanges")
:
Alternate ways of doing grglist(x)
and granges(x)
,
respectively.
In the code snippets below, x
is a GappedAlignmentPairs object.
show(x)
:
The number of lines displayed in the show method are controlled
with global options showHeadLines
and showTailLines
.
If the sum of the two options is greater than the object length
the full object is displayed. When these options are set they affect
the display of GRanges, GappedAlignments, GappedAlignmentPairs,
Ranges and XString objects.
H. Pages
ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools") galp <- readGappedAlignmentPairs(ex1_file, use.names=TRUE) galp length(galp) head(galp) head(names(galp)) first(galp) last(galp) last(galp, invert.strand=TRUE) left(galp) right(galp) seqnames(galp) strand(galp) head(ngap(galp)) table(isProperPair(galp)) seqlevels(galp) ## Rename the reference sequences: seqlevels(galp) <- sub("seq", "chr", seqlevels(galp)) seqlevels(galp) galp[[1]] unlist(galp) grglist(galp) # a GRangesList object grglist(galp, order.as.in.query=TRUE) stopifnot(identical(unname(elementLengths(grglist(galp))), ngap(galp) + 2L)) granges(galp) # a GRanges object introns(galp) # a GRangesList object stopifnot(identical(unname(elementLengths(introns(galp))), ngap(galp)))