To celebrate the release of GATK 4.0, we are giving away free credits for running the GATK4 Best Practices pipelines in FireCloud, our secure online analysis portal. It’s first come first serve, so sign up now to claim your free credits worth $250. Sponsored by Google Cloud. Learn more at https://software.broadinstitute.org/firecloud/documentation/freecredits

VCF file and allele frequency

Hello All,

I am using I am using GATK RNA-seq variant pipeline for finding muttaion/vatiants on the list of gene given in teh follwoing command line

java-1.7 -jar -Xincgc -Xmx1586M GenomeAnalysisTK-3.2-2.jar -T HaplotypeCaller --filter_reads_with_N_cigar -R human_genome37_gatk.fa -D dbsnp_137.hg19.vcf -I sample_split.bam -o sample.vcf -L mylist.intervals

And the resulting VCF files has for variants AF either 100 % or 50 % . It would be great if anyone would explain me what does AF means in INFO column from VCF file.
I AF means allele frequency or it has to be calculated separately from VCF file and if so how can I do it..???
Thank you

Tagged:

Best Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Accepted Answer

    The AF annotated by our tools is the theoretical allele frequency that corresponds to the genotype call made by the tool. If you want to calculate the empirical allele frequency you can use the values from the AD field. One way to do this is to extract the values using the VariantsToTable tool to produce a tab-delimited file, which you can then process with your preferred statistical tools.

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Accepted Answer

    The AF annotated by our tools is the theoretical allele frequency that corresponds to the genotype call made by the tool. If you want to calculate the empirical allele frequency you can use the values from the AD field. One way to do this is to extract the values using the VariantsToTable tool to produce a tab-delimited file, which you can then process with your preferred statistical tools.

  • AlvaAlva GermanyMember

    Ok, so i used the flowing command line ,

    java-1.7 -jar -Xincgc -Xmx1586M GenomeAnalysisTK-3.2-2.jar -R /human_genome37_gatk.fa -T VariantsToTable -V sample.vcf -F CHROM -F POS -F ID -F QUAL -F AD -o sample.table
    to retrieve the field.. But its throwing error as,

    ERROR
    ERROR MESSAGE: Missing field AD in vc variant at [VC variant @ chr1:14653 Q33.77 of type=SNP alleles=[C*, T] attr={AC=1, AF=0.500, AN=2, BaseQRankSum=0.720, ClippingRankSum=-1.380, DP=7, FS=0.000, MLEAC=1, MLEAF=0.500, MQ=60.00, MQ0=0, MQRankSum=-0.720, QD=4.82, ReadPosRankSum=1.380} GT=GT:AD:DP:GQ:PL 0/1:4,2:6:62:62,0,108

    But without AD as -F value its giving output table but I need AD..

    So one more question from there..AD is a comma separated value, as for example.. in the following variant,
    chr1 564598 . A G 14381.77 . AC=2;AF=1.00;AN=2;DP=342;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQ0=0;QD=29.56 GT:AD:DP:GQ:PL 1/1:0,342:342:99:14410,1029,0

    AD here is 0,342 : so the how is empirical Allele frequency calculated or which count is alternate varinat count and which is reference variant count.. .Could you please explain with this example...?
    Thank you

  • AlvaAlva GermanyMember

    Or putting the question in another way
    I ahve AD as GT:AD:DP:GQ:PL 1/1:2,366:368:99:14774,1067,0 in the column

    And for VAF teh formulate is variant reads / total reads and my question is which is variant read and which is ref read

  • AlvaAlva GermanyMember

    yeah I tried,
    with GF option and now I have the output file, but I have hard time understanding what is QUAL and AC columns means how one can infer Empirical Variant Allele frequency from these information .?
    the command line I used successfully,

    java-1.7 -jar -Xincgc -Xmx1586MGenomeAnalysisTK-3.2-2.jar -R human_genome37_gatk.fa -T VariantsToTable -V sample.vcf -F CHROM -F POS -F ID -F QUAL -GF -F AC -o sample.table
    and here is what i have as output
    head sample.table

    CHROM POS ID QUAL AC
    chr1 564598 rs6594028 13998.77 2
    chr1 564654 rs147404388 13921.77 2
    chr1 564862 rs1988726 646.77 2
    chr1 564868 rs1832728 646.77 2

    And from here how can I infer VAF (Variant allele frequency) I mean which values I need to use for that.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Are you looking for allele frequency in each sample, or overall in the population?

  • AlvaAlva GermanyMember
  • angezouangezou Member

    Hi @Sheila and @Geraldine_VdAuwera, i seem to have two separate DP valeues for one line, below is the output for one variant:

    1 680 rs794548358 C T 181.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-0.540;ClippingRankSum=-0.666;DB;DP=51;ExcessHet=3.0103;FS=10.533;MLEAC=1;MLEAF=0.500;MQ=44.74;MQRankSum=-1.972;QD=3.87;ReadPosRankSum=2.073;SOR=3.215 GT:AD:DP:GQ:PL 0/1:36,11:47:99:210,0,1046

    It says DP=51, but then later GT:AD:DP:GQ:PL 0/1:36,11:47:99:210,0,1046 lists DP=47, this only happens for some variants. Do you know why I am getting two different DP values and which one I should use? Also, why does allele frequency have to be AD/DP, if AD is 5, 35, can I also do 5/(35+5)? Thanks!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @angezou
    Hi,

    The difference in the two DPs has to do with filtering. Have a look at this page and this page for more information.

    The AD calculation is explained here.

    Also, one last article to help you! :smiley: https://www.broadinstitute.org/gatk/guide/article?id=6005

    I hope these help.

    -Sheila

  • WANGxiaojiWANGxiaoji ShanghaiMember

    @Geraldine_VdAuwera said:
    The AF annotated by our tools is the theoretical allele frequency that corresponds to the genotype call made by the tool. If you want to calculate the empirical allele frequency you can use the values from the AD field. One way to do this is to extract the values using the VariantsToTable tool to produce a tab-delimited file, which you can then process with your preferred statistical tools.

    Could you explain why the "theoretical allele frequency" is not equal to AC/AN?
    VCF reference document only says "use this when estimated from primary data, not called genotypes" to explain AF.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @WANGxiaoji
    Hi,

    Indeed the AF is equal to the AC/AN in GATK output. We differ a little from the VCF spec in that case.

    -Sheila

  • charlesbaudocharlesbaudo MissouriMember

    Hello,

    I'm attempting to estimate mitochondrial heteroplasmy using HaplotypeCaller on a single sample and I have a few questions.
    1) After perusing the forums it appears people have suggested using a ploidy of 100. Do you still believe this is a good recommendation?

    2) Can you please help me interpret the output to determine the alternate allele frequency for an individual site? I'm having difficulty following the above conversation and interpreting these results. I would like to make a statement similar to, "The alternate allele frequency for a given site ranges from 2-11%".

    Here is my command line:
    java -Xmx16g -jar GenomeAnalysisTK.jar -R /Volumes/Toshiba_NGS/Malassezia_sequences/TomDawson/MF1878.fasta -I /Volumes/Toshiba_NGS/Malassezia_sequences/TomDawson/1878final.bam -stand_call_conf 30 -stand_emit_conf 10 -ploidy 100 -o /Volumes/Toshiba_NGS/Malassezia_sequences/TomDawson/1878.raw.vcf -T HaplotypeCaller -minPruning 3 -nct 16

    and an output for two variants:

    1:
    mito 5092 . A T 35.17 my_snp_filter AC=1;AF=0.010;AN=100;BaseQRankSum=-0.187;ClippingRankSum=0.000;DP=1475;FS=0.000;MLEAC=1;MLEAF=0.010;MQ=60.00;MQRankSum=0.000;QD=0.02;ReadPosRankSum=-0.567;SOR=3.598 GT:AD:DP:GQ 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/1:1470,5:1475:53

    2:
    mito 27040 . T A 1561.98 my_snp_filter AC=8;AF=0.080;AN=100;BaseQRankSum=-9.864;ClippingRankSum=0.000;DP=1245;FS=93.623;MLEAC=8;MLEAF=0.080;MQ=44.95;MQRankSum=22.792;QD=1.52;ReadPosRankSum=-11.137;SOR=0.000 GT:AD:DP:GQ 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/1/1/1/1/1/1/1/1:929,98:1027:1

    Thanks,
    Charles

    Issue · Github
    by Sheila

    Issue Number
    1369
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
Sign In or Register to comment.