Use Htslib to extract specific region of vcf file

Question

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.?

Answer 1

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.