We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

SelectVariants on AF

estif74estif74 Saint Paul, MN, USAMember

Hi there,

I'm trying to do what I think should be pretty straightforward, which is to select sites in a VCF (generated by GATK) where the AF > 0.1. I have in my SelectVariants statement:

java -jar -Xmx8g /path/to/gatk-3.3-0/GenomeAnalysisTK.jar -T SelectVariants \
-R /path/to/ref/ref.fa \
-V /path/to/vcf/variants.vcf \
-select "AF > 0.1" \
-o /path/to/output.vcf

However I get the following message:

Invalid JEXL expression detected for select-0 with message ![0,8]: 'AF > 0.1;' > error

I'm wondering if the reason for this is because in my VCF there are some sites with multiple alleles (and hence multiple AFs) and is therefore an array? This is from my VCF header:

INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">

If this is the case, what would be the syntax to get GATK to look at the whole array of AFs? Or do I have to get rid of the multi-alleleic sites from my VCF? Or maybe I'm misinterpreting the error message?


Best Answer


  • estif74estif74 Saint Paul, MN, USAMember

    OK, well in the spirit of quickly answering my own question, and perhaps to benefit others, changing my select statement to:

    -select "AF.0 > 0.1"

    Seems to work great!

  • estif74estif74 Saint Paul, MN, USAMember

    Actually, I seem to still be having some trouble so it would be helpful if anyone has any ideas. It seems that while the GATK job will run, that it's not capturing all of the variants. Here is the new code:

    java -jar /path/to/GenomeAnalysisTK.jar -T SelectVariants \
    -R /path/to/ref.fa \
    -V /path/to/vcf/variants.vcf \
    -nt 8 \
    -env \
    -select "AF.0 > 0.1" \
    -o /path/to/output.vcf

    Here is an example of a variant that is in the input VCF file but doesn't make it into the output VCF file:

    chr8 26973544 . T C 22424.27 PASS AC=32;AF=0.571;AN=56;BaseQRankSum=-3.650e-01;ClippingRankSum=0.406;DP=640;FS=1.953;GQ_MEAN=184.95;GQ_STDDEV=171.88;InbreedingCoeff=0.0682;MQ=60.00;MQ0=0;MQRankSum=0.406;NCC=0;PG=4,0,1;POSITIVE_TRAIN_SITE;QD=22.16;ReadPosRankSum=0.023;SOR=0.861;VQSLOD=4.03;culprit=FS GT:AD:DP:FT:GQ:PL:PP 1/1:0,25:25:PASS:74:810,75,0:813,74,0 1/1:0,21:21:PASS:62:732,63,0:735,62,0 0/1:5,7:12:PASS:99:204,0,139:208,0,140 0/0:8,0:8:lowGQ:17:0,21,244:0,17,241 0/0:11,0:11:PASS:23:0,27,394:0,23,391 0/1:15,10:25:PASS:99:276,0,458:280,0,459 0/1:7,6:13:PASS:99:168,0,211:172,0,212 0/0:9,0:9:PASS:20:0,24,360:0,20,357 1/1:0,21:21:PASS:62:661,63,0:664,62,0 0/1:6,4:10:PASS:73:69,0,140:73,0,141 1/1:0,21:21:PASS:62:531,63,0:534,62,0 1/1:0,22:22:PASS:65:700,66,0:703,65,0 1/1:0,24:24:PASS:71:767,72,0:770,71,0 1/1:0,23:23:PASS:68:759,69,0:762,68,0 0/0:21,0:21:PASS:56:0,60,763:0,56,760 1/1:0,50:50:PASS:99:1742,150,0:1745,149,0 0/1:8,22:30:PASS:99:721,0,199:725,0,200 1/1:0,51:51:PASS:99:1797,153,0:1800,152,0 0/1:21,17:38:PASS:99:488,0,670:492,0,671 0/1:22,12:34:PASS:99:311,0,643:315,0,644 0/0:31,0:31:PASS:80:0,84,1260:0,80,1257 1/1:0,21:21:PASS:62:678,63,0:681,62,0 0/1:16,15:31:PASS:99:394,0,418:398,0,419 0/1:12,5:17:PASS:99:115,0,339:119,0,340 0/1:13,12:25:PASS:99:326,0,369:330,0,370 0/0:7,0:7:lowGQ:17:0,21,242:0,17,239 1/1:0,29:29:PASS:86:939,87,0:942,86,0 0/0:25,0:25:PASS:56:0,60,810:0,56,807

    Clearly the AF is > 0.1 so it should be included in the output VCF. Any thoughts?


  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    You were correct to suppose that the error you were getting was due to multi-allelic calls. Your workaround just accesses the first element in the array, so it's not doing the right thing for what you want, I think. Let me see if we have some workaround docs that I can dig up for you.

  • estif74estif74 Saint Paul, MN, USAMember

    Awesome thanks, any help you have would be much appreciated!

  • estif74estif74 Saint Paul, MN, USAMember

    Thanks for the info, never would have figured that out on my own by any stretch! It appears to be running smoothly and I will report back when the job is done. :smiley:


  • estif74estif74 Saint Paul, MN, USAMember

    Hey Geraldine that seemed to work great! Thanks so much for the help!!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Glad to hear it, @estif74 ! Credit for the magic formula goes to @gauthier, she'll be happy to hear that it has been helpful.

  • estif74estif74 Saint Paul, MN, USAMember

    sorry, meant to say @Geraldine_VdAuwera , @gauthier - you guys need a virtual tip jar :wink:

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭

    @estif74 said:
    you guys need a virtual tip jar :wink:

    Have you completely lost your marbles? I would be broke then :smiley: I just wanted to advertise for bcftools, which would be my weapon of choice for a task like this. I used to do this with awk :)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Well, Broad policy may not allow us to put a PayPal button on the forum homepage anyway, but we do appreciate thank you postcards!

Sign In or Register to comment.