Can some one help me in knowing algorithms for converting gVCF file for VCF file. As per https://sites.google.com/site/gvcftools/home/about-gvcf website, it can be done by removing non-vraiants regions from gVCF files but I am just wondering how to accomplish it means how to remove this ? Where to get coordinates for NON Variants region . I can write the scripts if once I know about how to accomplish it.
As of this post the newest version is v0.15. There's a tool called "extract_variants" in this package which can quickly convert the file to a standard variants-only VCF, for example:
You could whip up a script for either of these functions, but the gvcftools utilities can be faster/more convenient if you're dealing with a large number of samples. The code is also MIT licensed on github so you're welcome to fork and customize.
Hi Chris,
I was using extract_variants but its taking much time as compare to "grep -v 'BLOCKAVG' input.gvcf > output.vcf . So I am just wondering that if I may miss some rows if I used grep option or both are equally fine ?
Not all non-variant sites are converted to blocks -- you can have non-variants of length 1 which won't have a BLOCKAVG tag. If this matters to you, then the extract_variants option may be easier.
The only thing the grep command is doing is removing any lines that contain "BLOCKAVG". These lines are invariant regions. The remaining lines should be variant in at least one sample, otherwise they would be part of the block. The grep command is fast because it's not really doing much, where gvcftools is parsing each line in the file in a much more sophisticated way. I don't it's indicative of missing anything, but if you want to check the gvcftools output against grep, just run wc -l < file.vcf and make sure both files have the same number of lines.
Yes !! In fact I tried to use "extract_variants" but it's taking a large amount of time and in fact Its been 12 hour after running that and its still running. No doubt grep is much faster and easier also.
This doesn't sound right. The tool works as a filter so it might be waiting for stdin from you. If you're seeing it actively running for more than a few minutes (eg. by checking in top), then please report the issue if possible.
I tried running the file get_called_regions, I am getting nothing means 0kb file as ab output...I have checked it and repeat it many times still I ma not getting any output. Need some suggestions..and please add proper documentation also. Thanks
I tried your command but surprisingly I got the vcf file of same size and same number of lines means no filtration. When I check the .g.vcf file, I found that there is not keyword like 'BLOCKAVG' throughout the file. Why so? I generated the .g.vcf file in the same way as mentioned in gatk haplotype documentation. Need some suggestions. Thanks.
gVCF files are valid VCF files, so you should be able to use them in the same way as a normal VCF file.
Perhaps you could post an example (or link to) the files which you are using, so we can see in more detail exactly how they vary from standard VCF files.
You may get an example on the link provided in questions. gVCF files have additional information for non-vraiants regions also. I want to remove these information.
You could do it in Perl. Loop through line-by-line outputting only the lines where the alternative allele is not '.'. Probably not the most elegant solution, but it would work nevertheless.
Hi Chris, I was using extract_variants but its taking much time as compare to "grep -v 'BLOCKAVG' input.gvcf > output.vcf . So I am just wondering that if I may miss some rows if I used grep option or both are equally fine ?
Thanks Syed
Not all non-variant sites are converted to blocks -- you can have non-variants of length 1 which won't have a BLOCKAVG tag. If this matters to you, then the extract_variants option may be easier.
Thanks for responding to this. @syednajeebashraf: FYI Chris is the author of the gvcftools suite.
Thanks Chris, Ok !! I have already added grep version to my pipeline. But will try with extract_variants again. Will get back to you if need any help.
The only thing the grep command is doing is removing any lines that contain "BLOCKAVG". These lines are invariant regions. The remaining lines should be variant in at least one sample, otherwise they would be part of the block. The grep command is fast because it's not really doing much, where gvcftools is parsing each line in the file in a much more sophisticated way. I don't it's indicative of missing anything, but if you want to check the gvcftools output against grep, just run
wc -l < file.vcfand make sure both files have the same number of lines.Yes !! In fact I tried to use "extract_variants" but it's taking a large amount of time and in fact Its been 12 hour after running that and its still running. No doubt grep is much faster and easier also.
This doesn't sound right. The tool works as a filter so it might be waiting for stdin from you. If you're seeing it actively running for more than a few minutes (eg. by checking in top), then please report the issue if possible.
I tried running the file
get_called_regions, I am getting nothing means 0kb file as ab output...I have checked it and repeat it many times still I ma not getting any output. Need some suggestions..and please add proper documentation also. Thanks