GAlignmentPairs objects
The GAlignmentPairs class is a container for storing pairs of genomic alignments. These pairs are typically obtained by aligning paired-end reads to a reference genome or transcriptome.
A GAlignmentPairs object is a list-like object where each list element represents a pair of genomic alignment.
An alignment pair is made of a "first" and a "last"/"second" alignment, and is formally represented by a GAlignments object of length 2. In most applications, an alignment pair will represent an aligned paired-end read. In that case, the "first" member of the pair represents the alignment of the first end of the read (aka "first segment in the template", using SAM Spec terminology), and the "last" member of the pair represents the alignment of the second end of the read (aka "last segment in the template", using SAM Spec terminology).
In general, a GAlignmentPairs object will be created by loading
records from a BAM (or SAM) file containing aligned paired-end reads,
using the readGAlignmentPairs
function (see below).
Each element in the returned object will be obtained by pairing 2
records.
GAlignmentPairs(first, last, strandMode=1,
isProperPair=TRUE, names=NULL)
:
Low-level GAlignmentPairs constructor. Generally not used directly.
In the code snippets below, x
is a GAlignmentPairs object.
strandMode(x)
, strandMode(x) <- value
:
The strand mode is a per-object switch on GAlignmentPairs objects
that controls the behavior of the strand
getter.
More precisely, it indicates how the strand of a pair should be
inferred from the strand of the first and last alignments in the pair:
0: strand of the pair is always *.
1: strand of the pair is strand of its first alignment. This mode should be used when the paired-end data was generated using one of the following stranded protocols: Directional Illumina (Ligation), Standard SOLiD.
2: strand of the pair is strand of its last alignment. This mode should be used when the paired-end data was generated using one of the following stranded protocols: dUTP, NSR, NNSR, Illumina stranded TruSeq PE protocol.
These modes are equivalent to strandSpecific
equal 0, 1, and 2,
respectively, for the featureCounts
function defined in the
Rsubread package.
Note that, by default, the readGAlignmentPairs
function
sets the strand mode to 1 on the returned GAlignmentPairs object.
The function has a strandMode
argument to let the user set a
different strand mode. The strand mode can also be changed
any time with the strandMode
setter or with invertStrand
.
Also note that 3rd party programs TopHat2 and Cufflinks have a
--library-type
option to let the user specify which protocol
was used. Please refer to the documentation of these programs for more
information.
length(x)
:
Return the number of alignment pairs in x
.
names(x)
, names(x) <- value
:
Get or set the names on x
.
See readGAlignmentPairs
for how to automatically extract
and set the names when reading the alignments from a file.
first(x, real.strand=FALSE)
,
last(x, real.strand=FALSE)
,
second(x, real.strand=FALSE)
:
Get the "first" or "last"/"second" alignment for each alignment pair in
x
. The result is a GAlignments object of the same length
as x
.
If real.strand=TRUE
, then the strand is inverted on-the-fly
according to the strand mode currently set on the object (see
strandMode(x)
above). More precisely, if strandMode(x)
is 0, then the strand is set to * for the GAlignments object
returned by both, first()
and last()
.
If strandMode(x)
is 1, then the strand of the object returned
by last()
is inverted. If strandMode(x)
is 2, then the
strand of the object returned by first()
is inverted.
seqnames(x)
:
Get the sequence names of the pairs in x
i.e. the name of the
reference sequence for each alignment pair in x
. The sequence
name of a pair is the sequence name of the 2 alignments in the pair
if they are the same (concordant seqnames), or NA if they differ
(discordant seqnames).
The sequence names are returned in a factor-Rle
object that is parallel to x
, i.e. the i-th element in
the returned object is the sequence name of the i-th pair in x
.
strand(x)
:
Get the strand for each alignment pair in x
.
Obey strandMode(x)
above to infer the strand of a pair.
Return *
for pairs with discordant strand, or for
all pairs if strandMode(x)
is 0.
njunc(x)
:
Equivalent to njunc(first(x)) + njunc(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
GAlignmentPairs 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 GAlignmentPairs object.
x[i]
:
Return a new GAlignmentPairs object made of the selected
alignment pairs.
In the code snippets below, x
is a GAlignmentPairs object.
x[[i]]
:
Extract the i-th alignment pair as a GAlignments 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 GAlignments 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 GAlignmentPairs object.
granges(x, use.names=TRUE, use.mcols=FALSE,
on.discordant.seqnames=c("error", "drop", "split"))
,
ranges(x, use.names=TRUE, use.mcols=FALSE,
on.discordant.seqnames=c("error", "drop", "split"))
:
Return a GRanges object (for granges()
) or
IRanges) object (for ranges()
).
If x
contains no pairs with discordant seqnames, the
operation is guaranteed to be successful and to return an object
parallel to x
, that is, an object where the i-th element
is the range of the genomic region spanned by the i-th alignment in
x
(all gaps in the region are ignored).
If x
contains pairs with discordant seqnames, then an error
is raised, unless the on.discordant.seqnames
argument is set
to "drop"
or "split"
, in which case the pairs with
discordant seqnames are either dropped or represented with 2 genomic
ranges (or 2 ranges for ranges()
) in the returned object.
In that case, the returned object is NOT parallel to x
.
If use.names
is TRUE, then the names on x
(if any) are propagated to the returned object.
If use.mcols
is TRUE, then the metadata columns on x
(if any) are propagated to the returned object.
grglist(x, use.mcols=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
. The strand of
the returned ranges obeys the strand mode currently set on the
object (see strandMode(x)
above).
More precisely, if grl1
and grl2
are
grglist(first(x, real.strand=TRUE), order.as.in.query=TRUE)
and
grglist(last(x, real.strand=TRUE), order.as.in.query=TRUE)
,
respectively, then the i-th element in the returned GRangesList
object is c(grl1[[i]], grl2[[i]])
if strandMode(x)
is 0
or 1, or c(grl2[[i]], grl1[[i]])
if strandMode(x)
is 2.
Note that, if strandMode(x)
is 1 or 2, this results in the ranges
being in consistent order with the original "query template", that is,
being in the order defined by walking the "query template" from the
beginning to the end.
If use.names
is TRUE, then the names on x
(if any) are propagated to the returned object.
If use.mcols
is TRUE, then the metadata columns on x
(if any) are propagated to the returned object.
If drop.D.ranges
is TRUE, then deletions (Ds in the
CIGAR) are treated like junctions (Ns in the CIGAR), that is, the
ranges corresponding to deletions are dropped.
as(x, "GRanges")
, as(x, "IntegerRanges")
,
as(x, "GRangesList")
:
Alternate ways of doing
granges(x, use.names=TRUE, use.mcols=TRUE)
,
ranges(x, use.names=TRUE, use.mcols=TRUE)
, and
grglist(x, use.names=TRUE, use.mcols=TRUE)
, respectively.
as(x, "GAlignments")
:
Equivalent of unlist(x, use.names=TRUE)
.
In the code snippets below, x
is a GAlignmentPairs object.
show(x)
:
By default the show
method displays 5 head and 5 tail
elements. This can be changed by setting the global options
showHeadLines
and showTailLines
. If the object
length is less than (or equal to) the sum of these 2 options
plus 1, then the full object is displayed.
Note that these options also affect the display of GRanges
and GAlignments objects, as well as other Vector
derivatives defined in the IRanges and Biostrings packages
(e.g. IRanges and XStringSet objects).
Hervé Pagès
readGAlignmentPairs
for reading aligned paired-end
reads from a file (typically a BAM file) into a GAlignmentPairs
object.
GAlignments objects for handling aligned single-end reads.
makeGAlignmentPairs
for pairing the elements of a
GAlignments object into a GAlignmentPairs object.
junctions-methods for extracting and summarizing junctions from a GAlignmentPairs object.
coverage-methods for computing the coverage of a GAlignmentPairs object.
findOverlaps-methods for finding range overlaps between a GAlignmentPairs object and another range-based object.
seqinfo
in the GenomeInfoDb
package for getting/setting/modifying the sequence information
stored in an object.
The GRanges and GRangesList classes defined and documented in the GenomicRanges package.
library(Rsamtools) # for the ex1.bam file ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools") galp <- readGAlignmentPairs(ex1_file, use.names=TRUE, strandMode=1) galp length(galp) head(galp) head(names(galp)) first(galp) last(galp) # or second(galp) strandMode(galp) first(galp, real.strand=TRUE) last(galp, real.strand=TRUE) strand(galp) strandMode(galp) <- 2 first(galp, real.strand=TRUE) last(galp, real.strand=TRUE) strand(galp) seqnames(galp) head(njunc(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 strandMode(galp) <- 1 grglist(galp) ## Alternatively the strand mode can be changed with invertStrand(): invertStrand(galp) stopifnot(identical(unname(elementNROWS(grglist(galp))), njunc(galp) + 2L)) granges(galp) # a GRanges object
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.