<p>For iterating over reads in a region in a single BCF file, you want to use the iterator interface: <code>bcf_itr_queryi()</code> or <code>bcf_itr_querys()</code>, together with <code>bcf_itr_next()</code>. For a VCF file, it appears you need to use the similar TBX iterator interface: <code>tbx_itr_queryi()</code> or <code>tbx_itr_querys()</code>, together with <code>tbx_itr_next()</code>. See HTSlib's <em>tabix.c</em> for reasonably straightforward examples of both.</p>
<p>Alternatively you could use the routines in <em>htslib/synced_bcf_reader.h</em>, which encapsulates both and is what bcftools mostly uses.</p>

CC BY-SA 4.0
2021-05-03T19:30:48.487Z
2021-05-04T09:22:25.260Z
2021-05-04T09:22:25.260Z
John Marshall
<p>I am working in c++ using the htslib library to parse different vcf files as well as fasta files.</p>
<p>My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx_fetch_seq() and define a region of chr14:20000000-22000000.</p>
<p>Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.</p>
<pre><code>void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open(&quot;chr14.vcf&quot;, &quot;r&quot;);
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, &quot;20000000-22000000&quot;, 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout &lt;&lt; bcf_hdr_id2name(test_header, test_record-&gt;rid) &lt;&lt; &quot;\t&quot; &lt;&lt;
                          test_record-&gt;pos &lt;&lt; &quot;\t&quot; &lt;&lt;
                          test_record-&gt;n_allele &lt;&lt; &quot;\t&quot; &lt;&lt;
                          std::endl;
          }  
}
</code></pre>
<p>However, i doesn't work to extract all the variant sites in a specific region. I tried bcf_hdr_set_samples with the argument &quot;20000000-22000000&quot;, to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.</p>
<p>So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?</p>

<c++><htslib>
Use Htslib to extract specific region of vcf file
# Use Htslib to extract specific region of vcf file

(Asked by: username_0 on 2021-05-03)

I am working in c++ using the htslib library to parse different vcf files as well as fasta files.
My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx_fetch_seq() and define a region of chr14:20000000-22000000.
Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.
void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open("chr14.vcf", "r");
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, "20000000-22000000", 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout << bcf_hdr_id2name(test_header, test_record->rid) << "\t" <<
                          test_record->pos << "\t" <<
                          test_record->n_allele << "\t" <<
                          std::endl;
          }  
}

However, i doesn't work to extract all the variant sites in a specific region. I tried bcf_hdr_set_samples with the argument "20000000-22000000", to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.
So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?

**Tags:** c++, htslib

---

## 1 answers

### Answer by username_1 on 2021-05-03. Score: 4 | Accepted answer

For iterating over reads in a region in a single BCF file, you want to use the iterator interface: `bcf_itr_queryi()` or `bcf_itr_querys()`, together with `bcf_itr_next()`. For a VCF file, it appears you need to use the similar TBX iterator interface: `tbx_itr_queryi()` or `tbx_itr_querys()`, together with `tbx_itr_next()`. See HTSlib's *tabix.c* for reasonably straightforward examples of both.

Alternatively you could use the routines in *htslib/synced\_bcf\_reader.h*, which encapsulates both and is what bcftools mostly uses.


---

