Question: How to get the information of total number of homozygous and heterozygous non-reference from a multi sample VCF file?
0
gravatar for aneekbiotech
10 months ago by
India
aneekbiotech0 wrote:

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?

allele • 350 views
ADD COMMENTlink modified 10 months ago by Pierre Lindenbaum117k • written 10 months ago by aneekbiotech0

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.

ADD REPLYlink written 10 months ago by d-cameron2.0k

@ 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)..

ADD REPLYlink written 10 months ago by aneekbiotech0
0
gravatar for Pierre Lindenbaum
10 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum117k wrote:

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:

ADD COMMENTlink modified 10 months ago • written 10 months ago by Pierre Lindenbaum117k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1355 users visited in the last hour