Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?

Question

I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the following samtools command in C++:

samtools view name.bam -o name.sam scaffold:pos-pos

Since I must extract reads from thousands of regions, I do not want to iterate through the whole bam file each time I extract reads from a region. I am wondering if HTSlib can extract reads occurring in a specified scaffold and range without iterating through the whole bam. I found the function bam_itr_querys() in the htslib/sam.h documentation on GitHub, but don't know how to format the inputs for this function, or what the "idx" input is. I also saw somewhere that I would need to pair this function with bam_itr_next(), but don't know the format to do so.

EDIT: My bam files are sorted and indexed, if that matters.

Answer 1

#include <stdio.h>
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  hts_idx_t *idx = sam_index_load(fp, argv[1]);
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]);
  while (sam_itr_next(fp, itr, b) >= 0) {
    if (b->core.tid < 0) continue;
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos,
      b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)));
  }
  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
}

The usage would then be ./a.out alignments.bam chr1:123-456, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via sam_itr_querys() so it's not iterating over the entire file.

Answer 2

Have you considered using a sorted and indexed bam file? e.g.

samtools sort in.bam > out.sorted.bam
samtools index out.sorted.bam

This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the AlignmentFile.fetch() functionality.

Surely there is an analog in HTSlib API, though I don't know what to search to get it. It is in samtools as demonstrated here.

This answer seems to put the components together that you are interested in, among them sam_index_load(), sam_itr_querys().