text: user9118: 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.
  
  user37: #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.
  
original_text: 9118: <p>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++:</p>
  <pre><code>samtools view name.bam -o name.sam scaffold:pos-pos
  </code></pre>
  <p>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 &quot;idx&quot; 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.</p>
  <p>EDIT: My bam files are sorted and indexed, if that matters.</p>
   37: <pre class="lang-js prettyprint-override"><code>#include &lt;stdio.h&gt;
  #include &quot;htslib/sam.h&quot;
  int main(int argc, char *argv[])
  {
    samFile *fp = sam_open(argv[1], &quot;r&quot;);
    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 &lt; 3? &quot;.&quot; : argv[2]);
    while (sam_itr_next(fp, itr, b) &gt;= 0) {
      if (b-&gt;core.tid &lt; 0) continue;
      printf(&quot;%s\t%d\t%d\n&quot;, h-&gt;target_name[b-&gt;core.tid], b-&gt;core.pos,
        b-&gt;core.pos + bam_cigar2rlen(b-&gt;core.n_cigar, bam_get_cigar(b)));
    }
    bam_destroy1(b);
    hts_itr_destroy(itr);
    bam_hdr_destroy(h);
    sam_close(fp);
  }
  </code></pre>
  <p>The usage would then be <code>./a.out alignments.bam chr1:123-456</code>, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via <code>sam_itr_querys()</code> so it's not iterating over the entire file.</p>
  