How to get the information of total number of homozygous and heterozygous non-reference from a multi sample VCF file?

Question

[aneekbiotech]

How to get the information of total number of homozygous and heterozygous non-reference from a multi sample VCF file? As an example, I want the information like if a nonreference allele is B, then how many samples having AB and BB in that VCF file. Basically the number of samples. Is there any software tool available for that?

[d-cameron]

You can find that out by looking at the GT genotype field in the VCF file.

I recommend using the VariantAnnotation R package as you can do it in two lines: 1 to load the VCF, another to parse and summarise the VCF genotypes.

[aneekbiotech]

@ d-cameron,

Thank you very much for prompt reply, however I am not very familiar with the programming language. Is there any other way to do that, i.e. using available software like VCFtools, BCFtools, GATK. Also I want to add this information for each alternate allele as column in my annovar annotated variant file (in excel)..

Answer 1

[Pierre Lindenbaum]

using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

• stream all vcf
• convert to a stream of genotype
• filter out the HOM_REF
• convert to sampleName+"\t"+genotype-type
• count

ex:

java -jar dist/bioalcidaejdk.jar -e 'stream().flatMap(V->V.getGenotypes().stream()).filter(G->!G.isHomRef()).map(G->G.getSampleName()+"\t"+G.getType().name()).collect(Collectors.groupingBy(Function.identity(), Collectors.counting())).forEach((K,V)->println(K+"\t"+V));' src/test/resources/rotavirus_rf.vcf.gz

S3  HOM_VAR 8
S2  HOM_VAR 8
S1  HET 7
S2  HET 7
S3  HET 7
S5  HOM_VAR 8
S4  HET 7
S1  HOM_VAR 2
S4  HOM_VAR 7

---

another one: today, I wrote a GUI http://lindenb.github.io/jvarkit/VcfStatsJfx.html, one of the screens is a count of genotype-type per sample:

quickly wrote a java FX app displaying basic VCF stats : https://t.co/CxURRnKMre pic.twitter.com/chmS9ciuv9

— Pierre Lindenbaum (@yokofakun) April 9, 2018