GAlignments objects

Description

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

Details

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

Typically, a GAlignments 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 GAlignments 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 GAlignments 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 ?GAlignmentPairs for how to handle aligned paired-end reads).

Each element of a GAlignments object consists of:

• The name of the reference sequence. (This is the RNAME field in a SAM/BAM record.)

• The strand in the reference sequence to which the query is aligned. (This information is stored in the FLAG field in a SAM/BAM record.)

• The CIGAR string in the "Extended CIGAR format" (see the SAM Format Specifications for the details).

• The 1-based leftmost position/coordinate of the clipped query relative to the reference sequence. We will refer to it as the "start" of the query. (This is the POS field in a SAM/BAM record.)

• The 1-based rightmost position/coordinate of the clipped query relative to the reference sequence. We will refer to it as the "end" of the query. (This is NOT explicitly stored in a SAM/BAM record but can be inferred from the POS and CIGAR fields.) Note that all positions/coordinates are always relative to the first base at the 5' end of the plus strand of the reference sequence, even when the query is aligned to the minus strand.

• The genomic intervals between the "start" and "end" of the query that are "covered" by the alignment. Saying that the full [start,end] interval is covered is the same as saying that the alignment contains no junction (no N in the CIGAR). It is then considered to be a simple alignment. Note that a simple alignment can have mismatches or deletions (in the reference). In other words, a deletion (encoded with a D in the CIGAR) is NOT considered to introduce a gap in the coverage, but a junction is.

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

Optionally, a GAlignments 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:

• Access the information stored in a GAlignments object in a way that is independent from how the data are actually stored internally.

• How to create and manipulate a GAlignments object.

Constructor

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

Accessors

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

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

names(x), names(x) <- value: Get or set the names on x. See readGAlignments for how to automatically extract and set the names when reading the alignments from a file.

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).

njunc(x): Returns an integer vector of the same length as x containing the number of junctions (i.e. N operations in the CIGAR) in each alignment. Equivalent to unname(elementNROWS(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 GAlignments 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 GAlignments object.

granges(x, use.names=TRUE, use.mcols=FALSE), ranges(x, use.names=TRUE, use.mcols=FALSE):

Return a GRanges object (for granges()) or IRanges) object (for ranges()) parallel to x 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 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.names=TRUE, use.mcols=FALSE, order.as.in.query=FALSE, drop.D.ranges=FALSE), rglist(x, use.names=TRUE, use.mcols=FALSE, 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.

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.

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 junctions (N operations in the CIGAR), that is, the ranges corresponding to deletions are dropped.

See Details section above for more information.

as(x, "GRanges"), as(x, "Ranges"), as(x, "GRangesList"), as(x, "RangesList"): Alternate ways of doing granges(x, use.names=TRUE, use.mcols=TRUE), ranges(x, use.names=TRUE, use.mcols=TRUE), grglist(x, use.names=TRUE, use.mcols=TRUE), and rglist(x, use.names=TRUE, use.mcols=TRUE), respectively.

In the code snippet below, x is a GRanges object.

as(from, "GAlignments"): Creates a GAlignments object from a GRanges object. The metadata columns are propagated. cigar values are created from the sequence width unless a "cigar" metadata column already exists in from.

Subsetting and related operations

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

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

Combining

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

Other methods

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 GAlignmentPairs objects, as well as other objects defined in the IRanges and Biostrings packages (e.g. Ranges and DNAStringSet objects).

Author(s)

Herv<c3><83><c2><a9> Pag<c3><83><c2><a8>s and P. Aboyoun

References

http://samtools.sourceforge.net/

See Also

• readGAlignments for reading genomic alignments from a file (typically a BAM file) into a GAlignments object.

• GAlignmentPairs objects for handling aligned paired-end reads.

• junctions-methods for extracting and summarizing junctions from a GAlignments object.

• coverage-methods for computing the coverage of a GAlignments object.

• findOverlaps-methods for finding overlapping genomic alignments.

• 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.

• The CompressedIRangesList class defined and documented in the IRanges package.

Examples

library(Rsamtools)  # for the ex1.bam file
ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
gal <- readGAlignments(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(njunc(gal))
seqlevels(gal)

## Invert the strand:
invertStrand(gal)

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

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

## 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 junctions

## A confirmation that none of the alignments contains junctions (in
## other words, each alignment can be represented by a single genomic
## range on the reference):
stopifnot(all(njunc(gal) == 0))

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

## Unlike N operations, D operations don't introduce gaps:
ii <- grep("D", cigar(gal), fixed=TRUE)
gal[ii]
njunc(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]

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(GenomicAlignments)
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: GenomicRanges
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: Biostrings
Loading required package: XVector
Loading required package: Rsamtools
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/GenomicAlignments/GAlignments-class.Rd_%03d_medium.png", width=480, height=480)
> ### Name: GAlignments-class
> ### Title: GAlignments objects
> ### Aliases: class:GAlignments GAlignments-class GAlignments
> ###   updateObject,GAlignments-method length,GAlignments-method
> ###   names,GAlignments-method seqnames,GAlignments-method rname
> ###   rname,GAlignments-method strand,GAlignments-method
> ###   names<-,GAlignments-method seqnames<-,GAlignments-method rname<-
> ###   rname<-,GAlignments-method strand<-,GAlignments-method cigar
> ###   cigar,GAlignments-method qwidth qwidth,GAlignments-method
> ###   start,GAlignments-method end,GAlignments-method
> ###   width,GAlignments-method njunc njunc,GAlignments-method
> ###   seqinfo,GAlignments-method seqinfo<-,GAlignments-method
> ###   ranges,GAlignments-method granges,GAlignments-method
> ###   grglist,GAlignments-method rglist,GAlignments-method
> ###   coerce,GAlignments,Ranges-method coerce,GAlignments,GRanges-method
> ###   coerce,GAlignments,GRangesList-method
> ###   coerce,GAlignments,RangesList-method as.data.frame,GAlignments-method
> ###   coerce,GenomicRanges,GAlignments-method c,GAlignments-method
> ###   show,GAlignments-method
> ### Keywords: methods classes
> 
> ### ** Examples
> 
> library(Rsamtools)  # for the ex1.bam file
> ex1_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
> gal <- readGAlignments(ex1_file, param=ScanBamParam(what="flag"))
> gal
GAlignments object with 3271 alignments and 1 metadata column:
         seqnames strand       cigar    qwidth     start       end     width
            <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
     [1]     seq1      +         36M        36         1        36        36
     [2]     seq1      +         35M        35         3        37        35
     [3]     seq1      +         35M        35         5        39        35
     [4]     seq1      +         36M        36         6        41        36
     [5]     seq1      +         35M        35         9        43        35
     ...      ...    ...         ...       ...       ...       ...       ...
  [3267]     seq2      +         35M        35      1524      1558        35
  [3268]     seq2      +         35M        35      1524      1558        35
  [3269]     seq2      -         35M        35      1528      1562        35
  [3270]     seq2      -         35M        35      1532      1566        35
  [3271]     seq2      -         35M        35      1533      1567        35
             njunc |      flag
         <integer> | <integer>
     [1]         0 |        73
     [2]         0 |        73
     [3]         0 |       137
     [4]         0 |       137
     [5]         0 |       137
     ...       ... .       ...
  [3267]         0 |       137
  [3268]         0 |        73
  [3269]         0 |        83
  [3270]         0 |       147
  [3271]         0 |        83
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> ## ---------------------------------------------------------------------
> ## A. BASIC MANIPULATION
> ## ---------------------------------------------------------------------
> length(gal)
[1] 3271
> head(gal)
GAlignments object with 6 alignments and 1 metadata column:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     seq1      +         36M        36         1        36        36
  [2]     seq1      +         35M        35         3        37        35
  [3]     seq1      +         35M        35         5        39        35
  [4]     seq1      +         36M        36         6        41        36
  [5]     seq1      +         35M        35         9        43        35
  [6]     seq1      +         35M        35        13        47        35
          njunc |      flag
      <integer> | <integer>
  [1]         0 |        73
  [2]         0 |        73
  [3]         0 |       137
  [4]         0 |       137
  [5]         0 |       137
  [6]         0 |        73
  -------
  seqinfo: 2 sequences from an unspecified genome
> names(gal)  # no names by default
NULL
> seqnames(gal)
factor-Rle of length 3271 with 2 runs
  Lengths: 1482 1789
  Values : seq1 seq2
Levels(2): seq1 seq2
> strand(gal)
factor-Rle of length 3271 with 1394 runs
  Lengths:  8  1 63  2 12  1  2  2  6  1  4 ...  1  1  1  1  1 12  2  2  3  3
  Values :  +  -  +  -  +  -  +  -  +  -  + ...  +  -  +  -  +  -  +  -  +  -
Levels(3): + - *
> head(cigar(gal))
[1] "36M" "35M" "35M" "36M" "35M" "35M"
> head(qwidth(gal))
[1] 36 35 35 36 35 35
> table(qwidth(gal))

  33   34   35   36   40 
   6   37 2830  285  113 
> head(start(gal))
[1]  1  3  5  6  9 13
> head(end(gal))
[1] 36 37 39 41 43 47
> head(width(gal))
[1] 36 35 35 36 35 35
> head(njunc(gal))
[1] 0 0 0 0 0 0
> seqlevels(gal)
[1] "seq1" "seq2"
> 
> ## Invert the strand:
> invertStrand(gal)
GAlignments object with 3271 alignments and 1 metadata column:
         seqnames strand       cigar    qwidth     start       end     width
            <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
     [1]     seq1      -         36M        36         1        36        36
     [2]     seq1      -         35M        35         3        37        35
     [3]     seq1      -         35M        35         5        39        35
     [4]     seq1      -         36M        36         6        41        36
     [5]     seq1      -         35M        35         9        43        35
     ...      ...    ...         ...       ...       ...       ...       ...
  [3267]     seq2      -         35M        35      1524      1558        35
  [3268]     seq2      -         35M        35      1524      1558        35
  [3269]     seq2      +         35M        35      1528      1562        35
  [3270]     seq2      +         35M        35      1532      1566        35
  [3271]     seq2      +         35M        35      1533      1567        35
             njunc |      flag
         <integer> | <integer>
     [1]         0 |        73
     [2]         0 |        73
     [3]         0 |       137
     [4]         0 |       137
     [5]         0 |       137
     ...       ... .       ...
  [3267]         0 |       137
  [3268]         0 |        73
  [3269]         0 |        83
  [3270]         0 |       147
  [3271]         0 |        83
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> ## Rename the reference sequences:
> seqlevels(gal) <- sub("seq", "chr", seqlevels(gal))
> seqlevels(gal)
[1] "chr1" "chr2"
> 
> grglist(gal)  # a GRangesList object
GRangesList object of length 3271:
[[1]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1   [1, 36]      +

[[2]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames  ranges strand
  [1]     chr1 [3, 37]      +

[[3]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames  ranges strand
  [1]     chr1 [5, 39]      +

...
<3268 more elements>
-------
seqinfo: 2 sequences from an unspecified genome
> stopifnot(identical(unname(elementNROWS(grglist(gal))), njunc(gal) + 1L))
> granges(gal)  # a GRanges object
GRanges object with 3271 ranges and 0 metadata columns:
         seqnames       ranges strand
            <Rle>    <IRanges>  <Rle>
     [1]     chr1      [1, 36]      +
     [2]     chr1      [3, 37]      +
     [3]     chr1      [5, 39]      +
     [4]     chr1      [6, 41]      +
     [5]     chr1      [9, 43]      +
     ...      ...          ...    ...
  [3267]     chr2 [1524, 1558]      +
  [3268]     chr2 [1524, 1558]      +
  [3269]     chr2 [1528, 1562]      -
  [3270]     chr2 [1532, 1566]      -
  [3271]     chr2 [1533, 1567]      -
  -------
  seqinfo: 2 sequences from an unspecified genome
> rglist(gal)   # a CompressedIRangesList object
IRangesList of length 3271
[[1]]
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         1        36        36

[[2]]
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         3        37        35

[[3]]
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         5        39        35

...
<3268 more elements>
> stopifnot(identical(unname(elementNROWS(rglist(gal))), njunc(gal) + 1L))
> ranges(gal)   # an IRanges object
IRanges object with 3271 ranges and 0 metadata columns:
             start       end     width
         <integer> <integer> <integer>
     [1]         1        36        36
     [2]         3        37        35
     [3]         5        39        35
     [4]         6        41        36
     [5]         9        43        35
     ...       ...       ...       ...
  [3267]      1524      1558        35
  [3268]      1524      1558        35
  [3269]      1528      1562        35
  [3270]      1532      1566        35
  [3271]      1533      1567        35
> 
> ## Modify the number of lines in 'show'
> options(showHeadLines=3)
> options(showTailLines=2)
> gal
GAlignments object with 3271 alignments and 1 metadata column:
         seqnames strand       cigar    qwidth     start       end     width
            <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
     [1]     chr1      +         36M        36         1        36        36
     [2]     chr1      +         35M        35         3        37        35
     [3]     chr1      +         35M        35         5        39        35
     ...      ...    ...         ...       ...       ...       ...       ...
  [3270]     chr2      -         35M        35      1532      1566        35
  [3271]     chr2      -         35M        35      1533      1567        35
             njunc |      flag
         <integer> | <integer>
     [1]         0 |        73
     [2]         0 |        73
     [3]         0 |       137
     ...       ... .       ...
  [3270]         0 |       147
  [3271]         0 |        83
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> ## Revert to default
> options(showHeadLines=NULL)
> options(showTailLines=NULL)
> 
> ## ---------------------------------------------------------------------
> ## B. SUBSETTING
> ## ---------------------------------------------------------------------
> gal[strand(gal) == "-"]
GAlignments object with 1624 alignments and 1 metadata column:
         seqnames strand       cigar    qwidth     start       end     width
            <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
     [1]     chr1      -         35M        35        18        52        35
     [2]     chr1      -         35M        35       185       219        35
     [3]     chr1      -         34M        34       188       221        34
     [4]     chr1      -         35M        35       209       243        35
     [5]     chr1      -         35M        35       213       247        35
     ...      ...    ...         ...       ...       ...       ...       ...
  [1620]     chr2      -         33M        33      1523      1555        33
  [1621]     chr2      -         35M        35      1524      1558        35
  [1622]     chr2      -         35M        35      1528      1562        35
  [1623]     chr2      -         35M        35      1532      1566        35
  [1624]     chr2      -         35M        35      1533      1567        35
             njunc |      flag
         <integer> | <integer>
     [1]         0 |        89
     [2]         0 |       147
     [3]         0 |       121
     [4]         0 |        83
     [5]         0 |        83
     ...       ... .       ...
  [1620]         0 |        83
  [1621]         0 |       153
  [1622]         0 |        83
  [1623]         0 |       147
  [1624]         0 |        83
  -------
  seqinfo: 2 sequences from an unspecified genome
> gal[grep("I", cigar(gal), fixed=TRUE)]
GAlignments object with 27 alignments and 1 metadata column:
       seqnames strand       cigar    qwidth     start       end     width
          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
   [1]     chr1      +    18M5I12M        35       271       300        30
   [2]     chr1      +    14M5I17M        36       275       305        31
   [3]     chr1      +    11M5I19M        35       278       307        30
   [4]     chr2      +     9M2I24M        35       148       180        33
   [5]     chr2      +     9M2I24M        35       148       180        33
   ...      ...    ...         ...       ...       ...       ...       ...
  [23]     chr2      +    10M4I21M        35       775       805        31
  [24]     chr2      -     8M4I24M        36       777       808        32
  [25]     chr2      -     7M4I24M        35       778       808        31
  [26]     chr2      +     7M4I24M        35       778       808        31
  [27]     chr2      -     7M4I24M        35       778       808        31
           njunc |      flag
       <integer> | <integer>
   [1]         0 |        99
   [2]         0 |       163
   [3]         0 |       163
   [4]         0 |       163
   [5]         0 |       163
   ...       ... .       ...
  [23]         0 |       163
  [24]         0 |        83
  [25]         0 |        83
  [26]         0 |       163
  [27]         0 |       147
  -------
  seqinfo: 2 sequences from an unspecified genome
> gal[grep("N", cigar(gal), fixed=TRUE)]  # no junctions
GAlignments object with 0 alignments and 1 metadata column:
   seqnames strand       cigar    qwidth     start       end     width
      <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
       njunc |      flag
   <integer> | <integer>
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> ## A confirmation that none of the alignments contains junctions (in
> ## other words, each alignment can be represented by a single genomic
> ## range on the reference):
> stopifnot(all(njunc(gal) == 0))
> 
> ## Different ways to subset:
> gal[6]             # a GAlignments object of length 1
GAlignments object with 1 alignment and 1 metadata column:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chr1      +         35M        35        13        47        35
          njunc |      flag
      <integer> | <integer>
  [1]         0 |        73
  -------
  seqinfo: 2 sequences from an unspecified genome
> grglist(gal)[[6]]  # a GRanges object of length 1
GRanges object with 1 range and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1  [13, 47]      +
  -------
  seqinfo: 2 sequences from an unspecified genome
> rglist(gal)[[6]]   # a NormalIRanges object of length 1
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]        13        47        35
> 
> ## Unlike N operations, D operations don't introduce gaps:
> ii <- grep("D", cigar(gal), fixed=TRUE)
> gal[ii]
GAlignments object with 2 alignments and 1 metadata column:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chr2      -     9M1D26M        35       432       467        36
  [2]     chr2      -    13M1D22M        35      1501      1536        36
          njunc |      flag
      <integer> | <integer>
  [1]         0 |        83
  [2]         0 |       147
  -------
  seqinfo: 2 sequences from an unspecified genome
> njunc(gal[ii])
[1] 0 0
> grglist(gal[ii])
GRangesList object of length 2:
[[1]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames     ranges strand
         <Rle>  <IRanges>  <Rle>
  [1]     chr2 [432, 467]      -

[[2]] 
GRanges object with 1 range and 0 metadata columns:
      seqnames       ranges strand
  [1]     chr2 [1501, 1536]      -

-------
seqinfo: 2 sequences from an unspecified genome
> 
> ## qwidth() vs width():
> gal[qwidth(gal) != width(gal)]
GAlignments object with 29 alignments and 1 metadata column:
       seqnames strand       cigar    qwidth     start       end     width
          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
   [1]     chr1      +    18M5I12M        35       271       300        30
   [2]     chr1      +    14M5I17M        36       275       305        31
   [3]     chr1      +    11M5I19M        35       278       307        30
   [4]     chr2      +     9M2I24M        35       148       180        33
   [5]     chr2      +     9M2I24M        35       148       180        33
   ...      ...    ...         ...       ...       ...       ...       ...
  [25]     chr2      -     8M4I24M        36       777       808        32
  [26]     chr2      -     7M4I24M        35       778       808        31
  [27]     chr2      +     7M4I24M        35       778       808        31
  [28]     chr2      -     7M4I24M        35       778       808        31
  [29]     chr2      -    13M1D22M        35      1501      1536        36
           njunc |      flag
       <integer> | <integer>
   [1]         0 |        99
   [2]         0 |       163
   [3]         0 |       163
   [4]         0 |       163
   [5]         0 |       163
   ...       ... .       ...
  [25]         0 |        83
  [26]         0 |        83
  [27]         0 |       163
  [28]         0 |       147
  [29]         0 |       147
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> ## This MUST return an empty object:
> gal[cigar(gal) == "35M" & qwidth(gal) != 35]
GAlignments object with 0 alignments and 1 metadata column:
   seqnames strand       cigar    qwidth     start       end     width
      <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
       njunc |      flag
   <integer> | <integer>
  -------
  seqinfo: 2 sequences from an unspecified genome
> ## but this doesn't have too:
> gal[cigar(gal) != "35M" & qwidth(gal) == 35]
GAlignments object with 26 alignments and 1 metadata column:
       seqnames strand       cigar    qwidth     start       end     width
          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
   [1]     chr1      +    18M5I12M        35       271       300        30
   [2]     chr1      +    11M5I19M        35       278       307        30
   [3]     chr2      +     9M2I24M        35       148       180        33
   [4]     chr2      +     9M2I24M        35       148       180        33
   [5]     chr2      -     9M1D26M        35       432       467        36
   ...      ...    ...         ...       ...       ...       ...       ...
  [22]     chr2      +    10M4I21M        35       775       805        31
  [23]     chr2      -     7M4I24M        35       778       808        31
  [24]     chr2      +     7M4I24M        35       778       808        31
  [25]     chr2      -     7M4I24M        35       778       808        31
  [26]     chr2      -    13M1D22M        35      1501      1536        36
           njunc |      flag
       <integer> | <integer>
   [1]         0 |        99
   [2]         0 |       163
   [3]         0 |       163
   [4]         0 |       163
   [5]         0 |        83
   ...       ... .       ...
  [22]         0 |       163
  [23]         0 |        83
  [24]         0 |       163
  [25]         0 |       147
  [26]         0 |       147
  -------
  seqinfo: 2 sequences from an unspecified genome
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>