How to utilize the bai index file while parsing a BAM file

Question

[matejm469]

This is maybe a noobish question but I'm working on a program that visualizes some data from a SAM/BAM file, it also has an option to lookup certain reads by their names and stuff like that. I've been tasked to compare how much (if any) are some operations faster when performed on an indexed BAM file than on a simple SAM file but I don't know how would I go about utilizing a BAM index in my program. The program is written in java and I've been using htsjdk library for reading SAM/BAM files. I don't expect anyone to write my code for me, but some general pointers to how would I go about doing this would be much appreciated.

Thanks in advance

[rpolicastro]

I haven't used htsjdk, but their documentation shows support for BAM index files through their BAM reader API. Someone with more knowledge might come along and provide more info.

Answer 1

[Pierre Lindenbaum]

The bam index is just useful If you need to quickly get the reads in one or more defined region of the BAM. Eg:

How many reads are mapped on the gene 'X' on chr22

the file reader is quickly moved to the right section in the bam file without reading chr1, chr2, chr3... etc.. see https://javadoc.io/static/com.github.samtools/htsjdk/3.0.0/htsjdk/samtools/SamReader.html#query-htsjdk.samtools.QueryInterval:A-boolean-

For a queston like

How many reads is there in my bam file.

you don't need the BAI, just iterate over the whole bam. see: https://javadoc.io/static/com.github.samtools/htsjdk/3.0.0/htsjdk/samtools/SamReader.html#iterator--

To create a SamReader: https://javadoc.io/static/com.github.samtools/htsjdk/3.0.0/htsjdk/samtools/SamReaderFactory.html

edit: in htsjdk, you often don't need to specifify where is the bai file. It is automatically detected from the bam file.

[cmdcolin]

For "How many reads is there in my bam file." note that the BAI index does tell you this info, see samtools idxstats ("index stats") http://www.htslib.org/doc/samtools-idxstats.html

[Pierre Lindenbaum]

note that the BAI index does tell you this info, see samtools idxstats

how course it does :-) and samtools will always be faster for most simple jobs like this.

my example was only here to illustrate how to scan a bam using java.

[matejm469]

This might be a dumb question, but what exactly do you mean by "defined region of the BAM" ?

[Pierre Lindenbaum]

a genomic region defined by "chromosome:start-end"