Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-03T19:21:17.160Z
Answers.Text: My bam files are already sorted and indexed - I will check out this other answer
Answers.UserDisplayName: annabelperry
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-01T06:04:04.937Z
Answers.LastActivityDate: 2021-07-01T06:49:24.857Z
Answers.LastEditDate: 2021-07-01T06:49:24.857Z
Answers.UserDisplayName: Maximilian Press
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-03T19:19:37.643Z
Answers.Text: 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
Answers.UserDisplayName: annabelperry
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-03T20:19:25.387Z
Answers.Text: @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.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-03T23:18:06.480Z
Answers.Text: Thank you, @DevonRyan, that clears it up!
Answers.UserDisplayName: annabelperry
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-01T23:56:18.093Z
Answers.LastActivityDate: 2021-07-03T20:21:55.200Z
Answers.LastEditDate: 2021-07-03T20:21:55.200Z
Answers.UserDisplayName: user172818
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-06-30T20:14:41.490Z
Comments.Text: Have you looked at how samtools view works with a bed file of coordinates as input?
Comments.UserDisplayName: swbarnes2
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-07-01T19:40:48.243Z
Comments.Text: @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.
Comments.UserDisplayName: annabelperry
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-07-03T20:27:47.907Z
Comments.Text: @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.
Comments.UserDisplayName: Devon Ryan
Tags: <sam><samtools><c++><htslib>
Title: Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?
ThreadText: # 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**: @username_1 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**: @username_0 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
  
  ```js
  #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.
  
  > **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**: @username_0 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, @username_2, 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](https://www.biostars.org/p/151403/).
  
  [This answer](https://www.biostars.org/p/151053/) 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  
  
  ---
  
  