GappedAlignments-class {GenomicRanges} | R Documentation |
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.
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).
Most of the time, 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 ids (QNAME 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.
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".
The rest of this man page will focus on describing how to:
readGappedAlignments(file, format="BAM", ...)
:
Read a file as a GappedAlignments object.
The function is just a front-end that delegates to a format-specific
back-end function (any extra argument is 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 package
first).
In the code snippets below, x
is a GappedAlignments object.
length(x)
:
Returns the number of alignments in x
.
rname(x)
, rname(x) <- value
:
Gets or sets 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.
seqnames(x)
, seqnames(x) <- value
:
Same as rname(x)
and rname(x) <- value
.
strand(x)
, strand(x) <- value
:
Gets or sets 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).
grglist(x, drop.D.ranges = FALSE)
, granges(x)
,
rglist(x, drop.D.ranges = FALSE)
, ranges(x)
:
Returns either a GRangesList object, or a GRanges object,
or a RangesList object, or a Ranges
object of length length(x)
where each element represents the
regions in the reference to which a query is
aligned. If drop.D.ranges
is TRUE
for
either grglist
or rglist
, the ranges corresponding
to deletions in the CIGAR string are dropped, i.e., they are not
considered part of the alignment but are treated like the N (intron)
CIGAR element.
See Details section above for more information.
More precisely, the RangesList object returned by
rglist(x)
is a CompressedNormalIRangesList object,
and the Ranges object returned by ranges(x)
is
an IRanges object.
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 length length(x)
containing the number of gaps for each alignment.
Equivalent to elementLengths(rglist(x)) - 1L
.
seqinfo(x)
, seqinfo(x) <- value
:
Gets or sets the information about the underlying sequences.
value
must be a Seqinfo object.
seqlevels(x)
, seqlevels(x) <- value
:
Gets or sets the sequence levels.
seqlevels(x)
is equivalent to seqlevels(seqinfo(x))
or to levels(rname(x))
, those 2 expressions being
guaranteed to return indentical character vectors on a GappedAlignments
object. value
must be a character vector with no NAs.
seqlengths(x)
, seqlengths(x) <- value
:
Gets or sets 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
:
Gets or sets the circularity flags.
isCircular(x)
is equivalent to isCircular(seqinfo(x))
.
value
must be a named logical vector eventually with NAs.
In the code snippets below, x
is a GappedAlignments object.
x[i]
:
Returns a new GappedAlignments object made of the selected
alignments. i
can be a numeric or logical vector.
qnarrow(x, start=NA, end=NA, width=NA)
:
x
is a GappedAlignments object.
Returns 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).
narrow(x, start=NA, end=NA, width=NA)
:
x
is a GappedAlignments object.
Returns a new GappedAlignments object of the same length as x
describing the narrowed alignments. Unlike with qnarrow
now the start
/end
/width
arguments describe
the narrowing on the reference side, not the query side.
Like with qnarrow
, 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).
pintersect(x, y)
:
Either x
is a GappedAlignments object and y
is a
GRanges object or x
is a GRanges object and y
is a
GappedAlignments object.
Returns a new GappedAlignments object of the same length as the
GappedAlignments input arguments. Like with narrow
, the
resulting "parallel" intersection is with respect to the reference.
coverage(x, shift=0L, width=NULL, weight=1L)
:
Equivalent to coverage(as(x, "GRangesList"), shift, width, weight)
.
findOverlaps(query, subject, ...)
,
countOverlaps(query, subject, ...)
,
subsetByOverlaps(query, subject, ...)
,
match(x, table, nomatch=NA_integer_, incomparables=NULL)
,
x %in% table
:
query
or subject
or both are GappedAlignments objects.
findOverlaps(query, subject, ...)
is equivalent to
findOverlaps(grglist(query), subject, ...)
when
query
is a GappedAlignments object, or to
findOverlaps(query, grglist(subject), ...)
when
subject
is a GappedAlignments object, or to
findOverlaps(grglist(query), grglist(subject), ...)
when
both are GappedAlignments objects.
The same apply to countOverlaps(query, subject, ...)
and
subsetByOverlaps(query, subject, ...)
.
See ?`findOverlaps,GRangesList,GenomicRanges-method`
,
?`countOverlaps,GRangesList,GenomicRanges-method`
and
?`subsetByOverlaps,GRangesList,GenomicRanges-method`
for more information (in particular for descriptions of the
extra arguments and the returned object).
H. Pages and P. Aboyoun
http://samtools.sourceforge.net/
readBamGappedAlignments
,
GRangesList-class,
GRanges-class,
CompressedNormalIRangesList-class,
IRanges-class,
seqinfo
,
Seqinfo-class,
coverage
,
RleList-class,
pintersect,GRanges,GRanges-method,
findOverlaps,GRangesList,GenomicRanges-method,
countOverlaps,GRangesList,GenomicRanges-method,
subsetByOverlaps,GRangesList,GenomicRanges-method
library(Rsamtools) # the toy file we use below is in this package galn_file <- system.file("extdata", "ex1.bam", package="Rsamtools") galn <- readGappedAlignments(galn_file) galn ## --------------------------------------------------------------------- ## A. BASIC MANIPULATION ## --------------------------------------------------------------------- length(galn) head(galn) head(rname(galn)) seqlevels(galn) ## Rename the reference sequences: seqlevels(galn) <- sub("seq", "chr", seqlevels(galn)) seqlevels(galn) head(strand(galn)) head(cigar(galn)) head(qwidth(galn)) table(qwidth(galn)) grglist(galn) # a GRangesList object granges(galn) # a GRanges object rglist(galn) # a CompressedNormalIRangesList object ranges(galn) # an IRanges object stopifnot(identical(elementLengths(grglist(galn)), elementLengths(rglist(galn)))) head(start(galn)) head(end(galn)) head(width(galn)) head(ngap(galn)) ## --------------------------------------------------------------------- ## B. SUBSETTING ## --------------------------------------------------------------------- galn[strand(galn) == "-"] galn[grep("I", cigar(galn), fixed=TRUE)] galn[grep("N", cigar(galn), fixed=TRUE)] # no gaps ## A confirmation that all the queries map to the reference with no ## gaps: stopifnot(all(ngap(galn) == 0)) ## Different ways to subset: galn[6] # a GappedAlignments object of length 1 grglist(galn)[[6]] # a GRanges object of length 1 rglist(galn)[[6]] # a NormalIRanges object of length 1 ## Ds are NOT gaps: ii <- grep("D", cigar(galn), fixed=TRUE) galn[ii] ngap(galn[ii]) grglist(galn[ii]) ## qwidth() vs width(): galn[qwidth(galn) != width(galn)] ## This MUST return an empty object: galn[cigar(galn) == "35M" & qwidth(galn) != 35] ## but this doesn't have too: galn[cigar(galn) != "35M" & qwidth(galn) == 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(galn, 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(galn, start=21) stopifnot(identical(end(qnarrow(galn, start=21)), end(galn))) ## --------------------------------------------------------------------- ## D. coverage() ## --------------------------------------------------------------------- coverage(galn) ## --------------------------------------------------------------------- ## E. findOverlaps()/countOverlaps() ## --------------------------------------------------------------------- subject <- granges(galn)[1] subject ## Note the absence of query no. 9 (i.e. 'galn[9]') in this result: as.matrix(findOverlaps(galn, subject)) ## This is because findOverlaps()/countOverlaps() are strand specific: galn[8:10] countOverlaps(galn[8:10], subject) ## If this is not the desired behaviour, then you can set 'strand(galn)' ## to "*": strand(galn) <- "*" galn[8:10] countOverlaps(galn[8:10], subject) ## --------------------------------------------------------------------- ## F. ADVANCED OVERLAP EXAMPLES ## --------------------------------------------------------------------- subsetByOverlaps(galn, subject) table(match(galn, subject), useNA = "ifany") table(galn %in% subject)