This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Get Count Of 0|0 From A Gvcf

hi all, I have a gVCF file and i want to know which calls are 0|0 or wildtype homozygous. any ideas how this can be computed?

thanks.

vcf
• 4,105 views • link
updated by Chris Saunders Chris Saunders • written by win win

Re-call the variants jointly. If you do not have the BAM files you are out of luck. If you have the BAMs and don't want to re-call the variants you can check depth for each individual.

• 166 views • link
written by Zev.Kronenberg Zev.Kronenberg

2 answers

I don't have experience with this myself, but chaining commands such 'bcftools query' and 'awk' might do it:

bcftools query -f'%ALT\t%END\t%POS\n' <file> | awk '{if($1=="." && $2!=".")n+=($2-$3+1)}END{print n}'

The assumption here is that REF records in gVCF have ALT=. and the INFO/END annotation set. The command will print the total length of REF-only stretches. Depending on what your file looks like, it may need a bit more work.

Hope this helps

• 169 views • link
written by pd3 pd3

Using the latest gvcftools, you can run:

gzip -dc genome.vcf | ${GVCFTOOLS_BIN}/extract_variants --invert | ${GVCFTOOLS_BIN}/get_called_regions >| homref_regions.bed
• 0 views • link
written by Chris Saunders Chris Saunders

Log in to answer this question.

More posts like this