Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.
Attention:
We will be out of the office on November 11th and 13th 2019, due to the U.S. holiday(Veteran's day) and due to a team event(Nov 13th). We will return to monitoring the GATK forum on November 12th and 14th respectively. Thank you for your patience.

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?

Thanks,
Steve

Best Answer

Answers

  • 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?

    Thanks,
    Steve

  • 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:

    Best,
    Steve

  • 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.