<p>Have you considered using a sorted and indexed bam file? e.g.</p>
<pre><code>samtools sort in.bam &gt; out.sorted.bam
samtools index out.sorted.bam
</code></pre>
<p>This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the <code>AlignmentFile.fetch()</code> functionality.</p>
<p>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 <a href="https://www.biostars.org/p/151403/" rel="nofollow noreferrer">here</a>.</p>
<p><a href="https://www.biostars.org/p/151053/" rel="nofollow noreferrer">This answer</a> seems to put the components together that you are interested in, among them <code>sam_index_load()</code>, <code>sam_itr_querys()</code>.</p>

CC BY-SA 4.0
2021-07-03T19:21:17.160Z
My bam files are already sorted and indexed - I will check out this other answer
annabelperry
CC BY-SA 4.0
2021-07-01T06:04:04.937Z
2021-07-01T06:49:24.857Z
2021-07-01T06:49:24.857Z
Maximilian Press
<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>

CC BY-SA 4.0
2021-07-03T19:19:37.643Z
So I'm looking for a solution which doesn't require iterating through the whole bam file. Unless I am missing something, this solution requires that you iterate through the whole bam file and add "if-else" statements within the "for" loop to extract out your reads of interest. I am trying to avoid this, as it takes so much time. I'm looking for a solution which searches the binary bam file for reads in specific scaffolds/regions WITHOUT iterating through every read
annabelperry
CC BY-SA 4.0
2021-07-03T20:19:25.387Z
@annabelperry That isn't reading through the entire file. `sam_itr_querys` is accepting a string like `chr1:123-456` and then making an iterator over just those using the index (`idx`). This is all that samtools is doing internally. I've added a bit onto the end of this message to hopefully clarify that.
Devon Ryan
CC BY-SA 4.0
2021-07-03T23:18:06.480Z
Thank you, @DevonRyan, that clears it up!
annabelperry
CC BY-SA 4.0
2021-07-01T23:56:18.093Z
2021-07-03T20:21:55.200Z
2021-07-03T20:21:55.200Z
user172818
<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>

CC BY-SA 4.0
2021-06-30T20:14:41.490Z
Have you looked at how samtools view works with a bed file of coordinates as input?
swbarnes2
CC BY-SA 4.0
2021-07-01T19:40:48.243Z
@swbarnes2 I have not, but my issue is that I cannot use 'samtools view'. I should provide some context (I will edit the original question to include this info): in my first draft of my C++ program, I simply used the system() command to run 'samtools view' where needed. I am now wrapping the original program in an R package, and simply running the "system" command does not work. So, I need to run something with the same functionality as "samtools view" from within C++ itself.
annabelperry
CC BY-SA 4.0
2021-07-03T20:27:47.907Z
@annabelperry I can strongly recommend looking around at the internals of samtools if you need to work with htslib. Samtools is very well put together and I've found scrolling through its code to be a great way to gain insight into how to structure my own programs using htslib.
Devon Ryan
<sam><samtools><c++><htslib>
Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?
# Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?

(Asked by: username_0 on 2021-06-30)

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.

**Tags:** sam, samtools, c++, htslib
> **Comments:**
> **username_1**: Have you looked at how samtools view works with a bed file of coordinates as input?  
> **username_0**: @swbarnes2 I have not, but my issue is that I cannot use 'samtools view'. I should provide some context (I will edit the original question to include this info): in my first draft of my C++ program, I simply used the system() command to run 'samtools view' where needed. I am now wrapping the original program in an R package, and simply running the "system" command does not work. So, I need to run something with the same functionality as "samtools view" from within C++ itself.  
> **username_2**: @annabelperry I can strongly recommend looking around at the internals of samtools if you need to work with htslib. Samtools is very well put together and I've found scrolling through its code to be a great way to gain insight into how to structure my own programs using htslib.  

---

## 2 answers

### Answer by username_4 on 2021-07-01. Score: 3 | Accepted answer

#include 
#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.

> **Comments:**
> **username_0**: So I'm looking for a solution which doesn't require iterating through the whole bam file. Unless I am missing something, this solution requires that you iterate through the whole bam file and add "if-else" statements within the "for" loop to extract out your reads of interest. I am trying to avoid this, as it takes so much time. I'm looking for a solution which searches the binary bam file for reads in specific scaffolds/regions WITHOUT iterating through every read  
> **username_2**: @annabelperry That isn't reading through the entire file. `sam_itr_querys` is accepting a string like `chr1:123-456` and then making an iterator over just those using the index (`idx`). This is all that samtools is doing internally. I've added a bit onto the end of this message to hopefully clarify that.  
> **username_0**: Thank you, @DevonRyan, that clears it up!  

---

### Answer by username_3 on 2021-07-01. Score: 0

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().

> **Comments:**
> **username_0**: My bam files are already sorted and indexed - I will check out this other answer  

---

