This is a test version of Biostars. For the public version, visit https://www.biostars.org.
C++ how to read VCF file in indexed bgzip format?

How to quickly read VCF file in indexed bgzip format?

I am reading documents but get confused, there are libraries samtools, bcftools, htslib, vcftools... and I feel lost, which one should I use, could someone give me a quick C++ example, I just want to read line by line to get certain information of variants in a VCF file in indexed bgzip format.

samtools bcftools htslib vcftools bgzip
• 7,524 views • link
updated by rudy.d.arthur • rudy.d.arthur • written by Chen Sun Chen Sun

2 answers

it's clearly defined and readeable in the htslib source code: in the query_regions function

https://github.com/samtools/htslib/blob/master/tabix.c#L142

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

That is very helpful, is there any API documentation?

• 0 views • link
written by Chen Sun Chen Sun

Actually, if I just want to read bgzip VCF file line by line, do I really need index information? Thank you.

• 0 views • link
written by Chen Sun Chen Sun

no you can always 'just' use a gzip stream from gz lib : http://www.zlib.net/manual.html gz is compatible with the vcf/bgzf format

gzFile in= gzopen("input.vcf.gz", "rb");

it works much like stdio/ fopen

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

It seems there are differences between GBZF and GZIP file format. Should I use specific library or function to read BGZP file instead of this function for GZIP file?

• 0 views • link
written by Chen Sun Chen Sun

as I said bgzf is compatible with gzip. You can 'gzopen' a bgzf with gzlib but you cannot do a random-access of a bgzf without the htslib

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

To read bgzipped VCFs in C++ code I use htslib. Code to read a file looks something like

#include "htslib/hts.h"
extern "C" {
#include "htslib/synced_bcf_reader.h"
}
...
bcf_srs_t *sr =  bcf_sr_init() ;  //synced reader alloc
bcf_sr_add_reader (sr, input ); //open the file named 'input'
bcf1_t *line; //VCF line gets put in here
while(bcf_sr_next_line (sr)) { //loop through file
   line =  bcf_sr_get_line(sr, 0);  //read a line
   ...
}

The best thing to look at for API documentation is probably vcf.h

• 0 views • link
written by rudy.d.arthur • rudy.d.arthur

that is very helpful

• 0 views • link
written by Chen Sun Chen Sun

Log in to answer this question.

More posts like this