This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Hard time trying to calculate Allele Frequency and DP from Platypus

Hi I'm using Platypus for variant calling, but I'm having a hard time understanding the INFO header, most of which I cannot find in VCFv4.1 - ultimately what I'm trying to extract is DP and AF, however the header just does not make any sense to me nor can I find it in google. For example from the header can someone tell me how to extract DP and AD ( read count by allele)? 

BRF=0.0;FR=0.5000;HP=1;HapScore=1;MGOF=1;MMLQ=41;MQ=50.0;NF=19;NR=0;PP=675;QD=37.3563796077;SC=TTAAGAACACAGTGGCGCAGG;SbPval=1.0;Source=Platypus;TC=22;TCF=22;TCR=0;TR=19;WE=14472;WS=14454

GT:GL:GOF:GQ:NR:NV  0/1:-71.58,0.0,-5.98:1:60:22:19

thanks in advance!

rna-seq snv
• 5,000 views • link
updated by Kevin Blighe Kevin Blighe • written by simplitia • simplitia

Can you paste the VCF header? Remove the '##' if you can.

• 0 views • link
written by Kevin Blighe Kevin Blighe

@Kevin Blighe: yes thank you I'm still new to handling VCF files and your question was enough to answer my question since the header already contained the info I needed which was defining what each of those mean. For allele frequency I think this would be TR / TC since the header reads 

INFO=ID=TR,Number=.,Type=Integer,Description="Total number of reads containing this variant"
INFO=ID=NF,Number=.,Type=Integer,Description="Total number of forward reads containing this variant"
INFO=ID=TCF,Number=1,Type=Integer,Description="Total forward strand coverage at this locus"
INFO=ID=NR,Number=.,Type=Integer,Description="Total number of reverse reads containing this variant"
INFO=ID=TC,Number=1,Type=Integer,Description="Total coverage at this locus"
• 0 views • link
updated by Kevin Blighe Kevin Blighe • written by simplitia • simplitia

Yes, from that, I assume that the total read depth (DP) is equivalent to TC. I would then hope that TC were equivalent to NF+NR (?).

Finally, AF should be TR / TC, as you also mentioned.

• 0 views • link
written by Kevin Blighe Kevin Blighe

yes its good to check, I think you mean NF + NR = TR ; thanks again, super helpful.

• 0 views • link
written by simplitia • simplitia

1 answer

Further update:

You can add the AF, MAF, and other tags automatically via the BCFtools +fill-tags plugin.

See my answer, here, in particular part 4: A: How to use bcftools to calculate AF INFO field from AC and AN in VCF?

bcftools +fill-tags test.vcf
[W::bcf_hdr_check_sanity] PL should be declared as Number=G
5   135337248   .   CT  C   .   PASS    END=135337249;HOMLEN=3;HOMSEQ=TTT;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1  GT:AD   0/1:205,118
5   135337259   .   AG  A   .   PASS    END=135337260;HOMLEN=4;HOMSEQ=GGGG;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1 GT:AD   0/1:190,220
5   135337259   .   A   AG  .   PASS    END=135337259;HOMLEN=5;HOMSEQ=GGGGG;SVLEN=1;SVTYPE=INS;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1 GT:AD   0/1:192,71
5   135337264   .   GA  G   .   PASS    END=135337265;HOMLEN=1;HOMSEQ=A;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1    GT:AD   0/1:184,83
5   135337274   .   A   ATTATTGCATCAACTCCTCCGACATCTCTTCCCCTGCAAGAGTTCAGGCCCACAGGTTCTGGTGTGGGCTTGCTCAGCTGGAGGTAGCCTGAGGTGAGCTGGAG    .PASS   END=135337274;HOMLEN=23;HOMSEQ=TTATTGCATCAACTCCTCCGACA;SVLEN=103;SVTYPE=INS;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1    GT:AD   0/1:130,17

Kevin

• 1 views • link
written by Kevin Blighe Kevin Blighe

This is amazing, thank you :D

• 0 views • link
written by VBer VBer

Log in to answer this question.

More posts like this