AF in the final vcf

Hello,

I have my final vcf file after variant recalibration and I was wondering where the AF information comes from.
I suppose it comes from the training data. Right?

I am asking because I have too many variants with AF=1.00 or AF=0.00. and what I want is to extract the variants that have AF<0.15 according to the 1000 genome.

Can I do this through GATK? is the AF in my final file what I should be based on to do make the filtering?

Look forward to your reply.

Best,
Alexandra

Best Answers

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    The AF annotated by GATK is the theoretical allele frequency based on the genotype calls. It is determined by the variant caller and is not modified by the variant recalibration process.

    I'm not sure I understand what you're trying to do, can you give a detailed example?

  • alejandraalejandra spainMember

    Hello,

    Thank you for your response.

    If I understand correctly, the AF in the vcf file is the AF of the alternative allele.
    I want to select the variants that have a Minor Allele Frequency <0.20 or <0.15 according to the 1000 genome project. Can I do this through GATK? If not, do you have any idea of how I may do it?

    Also, I don't understand how the AF of some variants in the final file is 1.00 or 0.00. If this is correct, this is not a SNP. Right?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Can you show some example records from your file? It will be easier for us to explain what you're seeing if we can see what you're seeing.

    Also, please clarify what you mean by

    Minor Allele Frequency <0.20 or <0.15 according to the 1000 genome project.

    Do you mean you want to select variants from your dataset based on their MAF values in the 1000G project samples? Or based on their MAF in your own cohort?

  • alejandraalejandra spainMember

    Do you mean you want to select variants from your dataset based on their MAF values in the 1000G project samples? YES

    This is an example of my files. I have 29 samples.

    chr1 15250 . C CCAGGGT 284.54 . AC=2;AF=0.034;AN=58;BaseQRankSum=0.851;ClippingRankSum=1.66;DP=1380;FS=0.000;InbreedingCoeff=-0.0363;MLEAC=2;MLEAF=0.034;MQ=27.52;MQRankSum=1.42;QD=2.90;ReadPosRankSum=-1.955e+00;SOR=0.031 GT:AD:DP:GQ:PL 0/0:90,0:90:99:0,120,1800 0/0:68,0:68:99:0,102,1530 0/0:44,0:44:72:0,72,1080 0/0:57,0:57:99:0,114,1710 0/0:20,0:20:42:0,42,630 0/0:49,0:49:75:0,75,1125 0/0:67,0:67:99:0,108,1620 0/0:18,0:18:39:0,39,585 0/0:7,0:7:15:0,15,225 0/0:29,0:29:60:0,60,900 0/0:28,0:28:60:0,60,900 0/0:45,0:45:60:0,60,900 0/0:65,0:65:99:0,99,1485 0/0:30,0:30:48:0,48,720 0/0:44,0:44:78:0,78,1170 0/0:37,0:37:84:0,84,1260 0/0:54,0:54:93:0,93,1395 0/0:44,0:44:72:0,72,1080 0/0:24,0:24:48:0,48,720 0/0:33,0:33:66:0,66,990 0/0:46,0:46:99:0,99,1485 0/0:57,0:57:93:0,93,1395 0/0:16,0:16:36:0,36,540 0/1:34,9:43:99:275,0,1384 0/0:75,0:75:99:0,120,1800 0/1:50,5:55:59:59,0,2030 0/0:58,0:58:99:0,105,1575 0/0:74,0:74:99:0,114,1710 0/0:73,0:73:99:0,120,1800

    chr1 15274 rs201931625 A T 10123.84 PASS AC=56;AF=1.00;AN=56;DB;DP=360;FS=0.000;InbreedingCoeff=-0.0065;MLEAC=55;MLEAF=0.982;MQ=26.38;QD=29.60;SOR=10.341;VQSLOD=20.31;culprit=MQ GT:AD:DP:GQ:PL 1/1:0,29:29:87:836,87,0 1/1:0,8:8:24:244,24,0 1/1:0,5:8:15:147,15,0 1/1:0,17:17:50:505,50,0 1/1:0,4:4:12:114,12,0 1/1:0,8:8:24:236,24,0 1/1:0,27:29:81:781,81,0 1/1:0,4:4:12:115,12,0 ./.:0,0:0 1/1:0,6:6:18:177,18,0 1/1:0,6:6:18:182,18,0 1/1:0,4:4:12:114,12,0 1/1:0,14:14:42:395,42,0 1/1:0,6:6:18:179,18,0 1/1:0,10:10:30:294,30,0 1/1:0,9:9:27:266,27,0 1/1:0,11:11:33:325,33,0 1/1:0,8:8:24:239,24,0 1/1:0,5:5:15:147,15,0 1/1:0,5:5:15:154,15,0 1/1:0,14:14:42:416,42,0 1/1:0,15:15:45:432,45,0 1/1:0,5:5:15:158,15,0 1/1:0,15:15:45:426,45,0 1/1:0,30:33:90:888,90,0 1/1:0,18:20:54:532,54,0 1/1:0,19:19:59:580,59,0 1/1:0,20:20:60:596,60,0 1/1:0,20:20:60:597,60,0

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alejandra
    Hi,

    The AF field gives you the frequency that the alternate allele occurs in the genotypes of your samples. So, in the second example, all the sample genotypes are 1/1, giving an AF of 1.00.

    In your first example, two samples have 0/1 genotype. So, 2/58 = 0.034. The 58 comes from 29 samples * 2 possible alleles in the genotype.

    If you have two alternate alleles, the AF field will have two frequencies corresponding to the two alternate allele frequencies.

    If I understand correctly, you will use SelectVariants to select AF < "your threshold"

    This article may help as well.

    -Sheila

  • alejandraalejandra spainMember

    Hi Sheila,

    Thank you very much for your reply. Things are much clearer now.

    If I understand correctly, you will use SelectVariants to select AF < "your threshold" This article may help as well.

    I had already a look in the article before. This is not exactly what I want.
    I want to select the variants for which the minor allele frequency in the 1000 genome project is <0.20.
    I am looking for rare variants. Is there any way to do that through GATK. or do you have any idea of how I may extract those?

  • alejandraalejandra spainMember

    Thanks. I will try that.

    Also, I would like to ask 2 more questions that confused me.

    1. Considering the example that I sent. there is DP in the INFO field and also in the FORMAT field for each sample. What is the difference? and when I filter (SelectVariants) for DP >5, from which field is it taken? INFO or FORMAT?

    2. I have the idea that qual field >30 or >40 is a good variant but in some papers that use GATK protocol I saw that they use qual<30 to filter the variant. which one is the correct? or shall I filter for vqslod>0 that is the calibrated score. or both?

  • alejandraalejandra spainMember

    Hi, I noticed that the 1000G vcf is in assembly b37 and I have done my whole analysis with hg19. Do you think there is a way to avoid redoing the whole analysis using the b37. what do you suggest?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    If you can't find an hg19 version of the file you want to use, you can do a liftover to convert the file to the hg19 reference. See the Picard toolkit for the liftover tool.

  • alejandraalejandra spainMember

    Hi, can I use the 1000G_phase1.snps.high_confidence.hg19.sites.vcf from GATK bundle? it is from 1000 genome. Right?

  • alejandraalejandra spainMember

    HI Sheila,

    Thanks for your answer. related to that file 1000G.phase1.sites, I have the following question.

    I have some variants in my callset that do not exist in the 1000G site. One example is the chr1:14653:rs375086259:C:T in position 14653 that does exist in the NCBI database but does not exist in the 1000G site file from GATK bundle. More specifically, this file begins from chr1 position 51479.

    I noticed though that all the positions match to the file dbsnp138. I don't quite understand why the positions are so different since it is the same version of the genome.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator
    edited January 2016

    @alejandra
    Hi,

    I apologize, but I am not sure I understand your question. Are you asking why 1000Genomes VCF and dbSNP VCF have different sites of variation?

    -Sheila

  • alejandraalejandra spainMember

    Yes, it is the same version of the human genome reference. why do they have different sites?

    Also, for my whole analysis, I used the dbsnp.vcf file for the -known option. But I suppose I could use any vcf file that I considered fits better to the data. Is that right?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Why would you expect them to contain the same variant sites? Think about where they each come from (or read about it in the FAQs on resources if you don't know) and what they represent.

  • alejandraalejandra spainMember

    ok thanks. Now it makes sense. The 1000G in included in dbSNP but not the reverse.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
Sign In or Register to comment.