samtools: "sort by name" and fast random access

[cycomatto]

I understand that sorting a BAM file, e.g.

samtools sort my_file.bam my_file.sorted

will allow for faster random access. This makes sense, as viewing reads aligned to a certain region of the genome is much simpler if the reads are ordered by genome coordinates.

But if you sort a .bam file by read name, e.g.

samtools sort [B]-no[/B] my_file.bam my_file.sorted_by_name

then how does this allow faster random access?

I ask because there is no samtools functionality for observing reads by name - the "view" family of commands only accepts genome regions (coordinates).

How would you random access a .bam file by read name?

[maubp]

SAM/BAM files can be sorted by mapping position (normal), which allows efficient access to regions (e.g. chr1 from 1000 to 2000) via the BAI index, or sorted by read name (not commonly used).

If you want random access to reads by name, the BAM index file (BAI) is no help at all. You would need a separate index mapping read names to offsets. Sorting by read name would be helpful if you have multi-fragment reads (e.g. paired end reads), since then all the fragments would be together on disk so you'd only need one seek. If the reads were not sorted, you'd need to have a lookup table from read name to multiple offsets (one for each read fragment).

For BAM files you could use the same "virtual offset" approach as BAI uses (essentially a double offset, for a BGZF/GZIP block start, and the within block offset). You might find this blog post interesting about how indexing BGZF files (such as BAM files) can be done:

http://blastedbio.blogspot.com/2011/11/bgzf-blocked-bigger-better-gzip.html

and http://seqanswers.com/forums/showthread.php?t=15347

[cycomatto]

Sorting by read name would be helpful if you have multi-fragment reads (e.g. paired end reads), since then all the fragments would be together on disk so you'd only need one seek.

This is precisely my situation (given one side of a paired-end read, find its mate).


I am unclear why the following won't work: FIRST sorting by name, THEN index this name-sorted .bam.

Doesn't the "index" operation perform this offset for you? So if our .bam is already sorted by name, the the "index" operation will create the virtual offsets you wrote about?

[maubp]

Doesn't the "index" operation perform this offset for you? So if our .bam is already sorted by name, the the "index" operation will create the virtual offsets you wrote about?

The indexing provided by samtools/picard using the *.bai file ONLY works on BAM files sorted by mapping position.

I was talking about how as a programmer you can use your own alternative indexing, where you map read names to BAM/BGZF 'virtual offsets'.

[cycomatto]

So then, what is the point of sorting a .bam file by name? Why even include the option in samtools?

[kmcarr]

So then, what is the point of sorting a .bam file by name? Why even include the option in samtools?

If you want to process the entire file sequentially and you need to operate on read pairs as a single unit having reads sorted by name (and hence read pairs will be adjacent in the file) greatly reduces complexity. An example of a program which works this way is htseq-count.

[ajay87]

htseq-count also requires the file to be sorted

[xied75]

So then, what is the point of sorting a .bam file by name? Why even include the option in samtools?

samtools fixmate