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