Rsamtools

The BAM / SAM file format

The SAM format is a text based representation of alignments. The BAM format is a binary version of SAM which is smaller and much faster. In general, always work with BAM.

The format is quite complicated, which means the R representation is also a bit complicated. This complication happens because of the following features of the file format

• It may contain unaligned sequences.
• Each sequence may be aligned to multiple locations.
• It supports spliced (split) alignments.
• It may contain reads from multiple samples.

We will not attempt to fully understanding all the intricacies of this format.

A BAM file can be sorted in multiple ways. If it is sorted according to genomic location and if it is also “indexed” it is possible to retrieve all reads mapping to a genomic location, something which can be very useful. In Rsamtools this is done by the functions sortBam() and indexBam().

Reading in parts of the file

BAM files can be extremely big and it is there often necessary to read in parts of the file. You can do this in different ways

1. Read a set number of records (alignments).
2. Only read alignments satisfying certain criteria.
3. Only read alignments in certain genomic regions.

Let us start with the first of this. By specifying yieldSize when you use BamFile(), every invocation of scanBam() will only read yieldSize number of alignments. You can then invoke scanBam() again to get the next set of alignments; this requires you to open() the file first (otherwise you will keep read the same alignments).

yieldSize(bamFile) <- 1
open(bamFile)
scanBam(bamFile)[[1]]$seq

##   A DNAStringSet instance of length 1
##     width seq
## [1]    36 CACTAGTGGCTCATTGTAAATGTGTGGTTTAACTCG

scanBam(bamFile)[[1]]$seq

##   A DNAStringSet instance of length 1
##     width seq
## [1]    35 CTAGTGGCTCATTGTAAATGTGTGGTTTAACTCGT

## Cleanup
close(bamFile)
yieldSize(bamFile) <- NA

The other two ways of reading in parts of a BAM file is to use ScanBamParams(), specifically the what and which arguments.

gr <- GRanges(seqnames = "seq2",
              ranges = IRanges(start = c(100, 1000), end = c(1500,2000)))
params <- ScanBamParam(which = gr, what = scanBamWhat())
aln <- scanBam(bamFile, param = params)
names(aln)

## [1] "seq2:100-1500"  "seq2:1000-2000"

head(aln[[1]]$pos)

## [1] 66 68 68 72 73 77

Notice how the pos is less than what is specified in the which argument; this is because the alignments overlap the which argument. The what=scanBamWhat() tells the function to read everything. Often, you may not be interested in the actual sequence of the read or its quality scores. These takes up a lot of space so you may consider disabling reading this information.

BAM summary

Sometimes you want a quick summary of the alignments in a BAM file:

quickBamFlagSummary(bamFile)

##                                 group |    nb of |    nb of | mean / max
##                                    of |  records |   unique | records per
##                               records | in group |   QNAMEs | unique QNAME
## All records........................ A |     3307 |     1699 | 1.95 / 2
##   o template has single segment.... S |        0 |        0 |   NA / NA
##   o template has multiple segments. M |     3307 |     1699 | 1.95 / 2
##       - first segment.............. F |     1654 |     1654 |    1 / 1
##       - last segment............... L |     1653 |     1653 |    1 / 1
##       - other segment.............. O |        0 |        0 |   NA / NA
## 
## Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.
## Indentation reflects this.
## 
## Details for group M:
##   o record is mapped.............. M1 |     3271 |     1699 | 1.93 / 2
##       - primary alignment......... M2 |     3271 |     1699 | 1.93 / 2
##       - secondary alignment....... M3 |        0 |        0 |   NA / NA
##   o record is unmapped............ M4 |       36 |       36 |    1 / 1
## 
## Details for group F:
##   o record is mapped.............. F1 |     1641 |     1641 |    1 / 1
##       - primary alignment......... F2 |     1641 |     1641 |    1 / 1
##       - secondary alignment....... F3 |        0 |        0 |   NA / NA
##   o record is unmapped............ F4 |       13 |       13 |    1 / 1
## 
## Details for group L:
##   o record is mapped.............. L1 |     1630 |     1630 |    1 / 1
##       - primary alignment......... L2 |     1630 |     1630 |    1 / 1
##       - secondary alignment....... L3 |        0 |        0 |   NA / NA
##   o record is unmapped............ L4 |       23 |       23 |    1 / 1

BamViews

Instead of reading a single file, it is possible to construct something called a BamViews, a link to multiple files. In many ways, it has the same Views functionality as other views. A quick example should suffice, first we read everything;

bamView <- BamViews(bamPath)
aln <- scanBam(bamView)
names(aln)

## [1] "ex1.bam"

This gives us an extra list level on the return object; first level is files, second level is ranges.

We can also set bamRanges() on the BamViews to specify that only certain ranges are read; this is similar to setting a which argument to ScanBamParams().

bamRanges(bamView) <- gr
aln <- scanBam(bamView)
names(aln)

## [1] "ex1.bam"

names(aln[[1]])

## [1] "seq2:100-1500"  "seq2:1000-2000"

countBam

Sometimes, all you want to do is count… use countBam() instead of scanBam().