How To Get Count Of 0|0 From A Gvcf

Question

[win]

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.

[Zev.Kronenberg]

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.

Answer 1

[pd3]

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

Answer 2

[Chris Saunders]

Using the latest gvcftools, you can run:

gzip -dc genome.vcf | ${GVCFTOOLS_BIN}/extract_variants --invert | ${GVCFTOOLS_BIN}/get_called_regions >| homref_regions.bed