Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
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 for a Broad Institute event from Dec 10th to Dec 11th 2019. We will be back to monitor the GATK forum on Dec 12th 2019. In the meantime we encourage you to help out other community members with their queries.
Thank you for your patience!

Variantfiltration not recognizing first position annotations

Dear GATK staffs,

I tried to do Variantfiltration of my result from GenotypeGVCFs so that i can combined the variant into a single vcf file before pipe into VQSR. This is because my 31 vcf files have very few variant, ranging from 100 to 3000 which makes the Gausian model less reliable(full of FP showed in the graph as well). Feel free to advice me if i am doing silly mistakes in this workflow.

However, in both enxountering of Variantfiltration (using my result from GVCFs and result after VQSR) , i faced the same problem.

17:01:16.119 WARN JexlEngine - ![0,14]: 'ReadPosRankSum < -8.0 || MQRankSum < -12.5 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0;' undefined variable ReadPosRankSum
17:01:16.121 WARN JexlEngine - ![0,14]: 'ReadPosRankSum < -8.0 || MQRankSum < -12.5 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0;' undefined variable ReadPosRankSum
17:01:16.125 WARN JexlEngine - ![0,14]: 'ReadPosRankSum < -8.0 || MQRankSum < -12.5 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0;' undefined variable ReadPosRankSum
17:01:16.126 WARN JexlEngine - ![0,14]: 'ReadPosRankSum < -8.0 || MQRankSum < -12.5 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0;' undefined variable ReadPosRankSum

When i swapped ReadposRankSum with other annotation, those that placed in the first position in the beginning will be reported undefined. Hope that i can find an answer here.

My commands:
for file in *.vcf.gz; do gatk VariantFiltration -R $reference_dir -O ${file%%.vcf.gz}_filtered.vcf.gz -V $file --filter-name "snps_filter" --filter-expression "ReadPosRankSum < -8.0 || MQRankSum < -12.5 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0" ; done

Answers

  • wlaiwlai Member
    After several attempts such as running the the annotations with parenthesis and independenly, i think variantfiltration doesnt recognised annotations under INFO column .
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @wlai

    From the error it looks like the variants do not have ReadPosRankSum annotation. Can you please post a few records of the vcf file you are trying to filter.

  • wlaiwlai Member
    hi @bhanuGandham ,
    thanks for the reply. This is the current file left with me( after VQSR). Previous file has been removed but i think they showed the similar headers.

    #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT LB191
    chr1 169464512 . G T 26.63 VQSRTrancheSNP55.90to60.00+ AC=0;AF=0.00;AN=0;BaseQRankSum=0.00;DP=0;Exce
    ssHet=3.0103;FS=0.000;InbreedingCoeff=-0.1559;MQ=60.00;MQRankSum=0.00;QD=3.80;ReadPosRankSum=-3.660e-01;SOR=0.446;VQSLOD=15.31;culpri
    t=QD GT:AD:DP:PL ./.:0,0:0:0,0,0
    chr1 169464882 . T C 25.27 VQSRTrancheSNP55.90to60.00+ AC=0;AF=0.00;AN=2;BaseQRankSum=0.00;DP=2;Exce
    ssHet=3.0103;FS=0.000;InbreedingCoeff=-0.0992;MQ=60.00;MQRankSum=0.00;POSITIVE_TRAIN_SITE;QD=4.21;ReadPosRankSum=-1.834e+00;SOR=0.693
    ;VQSLOD=16.67;culprit=QD GT:AD:DP:GQ:PL 0/0:2,0:2:6:0,6,39
    chr1 169465177 . G T 21.76 VQSRTrancheSNP55.90to60.00+ AC=0;AF=0.00;AN=0;BaseQRankSum=0.00;DP=0;Exce
    ssHet=3.0103;FS=0.000;InbreedingCoeff=-0.1681;MQ=60.00;MQRankSum=0.00;QD=2.72;ReadPosRankSum=0.272;SOR=0.693;VQSLOD=15.30;culprit=QD
    GT:AD:DP:PL ./.:0,0:0:0,0,0


    I have some doubts with my result from ApplyVQSR, which leads to this post as i thought this could be the cause. However, i failed to post them for unknown reasons(asking me to be here a little longer) for more than 3 weeks but it seems fine to me to post other posts. Any suggestion to resolve this?
    Hope it is okay for me to ask you my ultimate question at here. I have stucked at this point for a long while....

    In total I have 45108 variants from 31 exome vcf files. After VariantRecalibrator, my tranche specificity shows the best Ti/Tv ratio at 57.8 which is very different from several examples from the tutorial(99.9). Is it normal and should i proceed with ApplyVQSR with -tranche cutoff level at 55.7? I have read several tutorials and video which tells me it is alright to lower down the value but i have not yet seen anyone has lowered it to below 80.

    If this is not normal, should i switch to hardfilter?

    Command:

    gatk --java-options "-Xmx48g -Xms48g" VariantRecalibrator -R $reference_dir -V tmp_sitesonly.vcf.gz -O tmp_sitesonly_recal_snps.vcf.gz --tranches-file tmp_sitesonly.snps.tranches -an MQ -an DP -an QD -an FS -an ReadPosRankSum --max-gaussians 4 --trust-all-polymorphic -tranche 60 -tranche 55.9 -tranche 55.8 -tranche 55.7 -tranche 55.5 -tranche 55 -mode SNP --rscript-file tmp_Merged_snps.plot.R --resource hapmap,known=false,training=true,truth=true,prior=15.0:$hapmap --resource omni,known=false,training=true,truth=false,prior=12.0:$omni --resource 1000G,known=false,training=true,truth=false,prior=10.0:$oneksnp --resource dbsnp,known=true,training=false,truth=false,prior=2.0:$dbsnp

    Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xmx48g -Xms48g -jar /home/cruadmin01/installations/gatk-4.0.12.0/gatk-package-4.0.12.0-local.jar VariantRecalibrator -R /home/Graceca/references/hs38/hs38DH.fa -V tmp_sitesonly.vcf.gz -O tmp_sitesonly_recal_snps.vcf.gz --tranches-file tmp_sitesonly.snps.tranches -an MQ -an DP -an QD -an FS -an ReadPosRankSum --max-gaussians 4 --trust-all-polymorphic -tranche 60 -tranche 55.9 -tranche 55.8 -tranche 55.7 -tranche 55.5 -tranche 55 -mode SNP --rscript-file tmp_Merged_snps.plot.R --resource hapmap,known=false,training=true,truth=true,prior=15.0:/home/Graceca/references/GATK/hapmap_3.3.hg38.vcf.gz --resource omni,known=false,training=true,truth=false,prior=12.0:/home/Graceca/references/GATK/1000G_omni2.5.hg38.vcf.gz --resource 1000G,known=false,training=true,truth=false,prior=10.0:/home/Graceca/references/GATK/1000G_phase1.snps.high_confidence.hg38.vcf.gz --resource dbsnp,known=true,training=false,truth=false,prior=2.0:/home/Graceca/references/GATK/dbsnp_146.hg38.vcf.gz
    16:13:44.468 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/cruadmin01/installations/gatk-4.0.12.0/gatk-package-4.0.12.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
    16:13:46.193 INFO VariantRecalibrator - ------------------------------------------------------------
    16:13:46.193 INFO VariantRecalibrator - The Genome Analysis Toolkit (GATK) v4.0.12.0
    16:13:46.194 INFO VariantRecalibrator - For support and documentation go to https://software.broadinstitute.org/gatk/
    16:13:46.194 INFO VariantRecalibrator - Executing as [email protected] on Linux v3.13.0-103-generic amd64
    16:13:46.194 INFO VariantRecalibrator - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_101-b13
    16:13:46.194 INFO VariantRecalibrator - Start Date/Time: 22 March, 2019 4:13:44 PM SGT
    16:13:46.194 INFO VariantRecalibrator - ------------------------------------------------------------
    16:13:46.195 INFO VariantRecalibrator - ------------------------------------------------------------
    16:13:46.195 INFO VariantRecalibrator - HTSJDK Version: 2.18.1
    16:13:46.195 INFO VariantRecalibrator - Picard Version: 2.18.16
    16:13:46.195 INFO VariantRecalibrator - HTSJDK Defaults.COMPRESSION_LEVEL : 2
    16:13:46.195 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
    16:13:46.196 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
    16:13:46.196 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
    16:13:46.196 INFO VariantRecalibrator - Deflater: IntelDeflater
    16:13:46.196 INFO VariantRecalibrator - Inflater: IntelInflater
    16:13:46.196 INFO VariantRecalibrator - GCS max retries/reopens: 20
    16:13:46.196 INFO VariantRecalibrator - Requester pays: disabled
    16:13:46.196 INFO VariantRecalibrator - Initializing engine
    16:13:46.593 INFO FeatureManager - Using codec VCFCodec to read file file:///home/Graceca/references/GATK/hapmap_3.3.hg38.vcf.gz
    16:13:46.733 INFO FeatureManager - Using codec VCFCodec to read file file:///home/Graceca/references/GATK/1000G_omni2.5.hg38.vcf.gz
    16:13:46.833 INFO FeatureManager - Using codec VCFCodec to read file file:///home/Graceca/references/GATK/1000G_phase1.snps.high_confidence.hg38.vcf.gz
    16:13:46.903 INFO FeatureManager - Using codec VCFCodec to read file file:///home/Graceca/references/GATK/dbsnp_146.hg38.vcf.gz
    16:13:46.957 INFO FeatureManager - Using codec VCFCodec to read file file:///mnt1/fastq/hs38dr/tmp_sitesonly.vcf.gz
    16:13:47.041 WARN IndexUtils - Feature file "/home/Graceca/references/GATK/dbsnp_146.hg38.vcf.gz" appears to contain no sequence dictionary. Attempting to retrieve a sequence dictionary from the associated index file
    16:13:47.148 INFO VariantRecalibrator - Done initializing engine
    16:13:47.150 INFO TrainingSet - Found hapmap track: Known = false Training = true Truth = true Prior = Q15.0
    16:13:47.150 INFO TrainingSet - Found omni track: Known = false Training = true Truth = false Prior = Q12.0
    16:13:47.150 INFO TrainingSet - Found 1000G track: Known = false Training = true Truth = false Prior = Q10.0
    16:13:47.150 INFO TrainingSet - Found dbsnp track: Known = true Training = false Truth = false Prior = Q2.0
    16:13:47.206 INFO ProgressMeter - Starting traversal
    16:13:47.206 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
    16:13:52.601 INFO ProgressMeter - chr20:44357841 0.1 49559 551268.1
    16:13:52.601 INFO ProgressMeter - Traversal complete. Processed 49559 total variants in 0.1 minutes.
    16:13:52.616 INFO VariantDataManager - MQ: mean = 59.81 standard deviation = 1.40
    16:13:52.629 INFO VariantDataManager - DP: mean = 90.09 standard deviation = 369.36
    16:13:52.635 INFO VariantDataManager - QD: mean = 23.68 standard deviation = 7.09
    16:13:52.640 INFO VariantDataManager - FS: mean = 0.16 standard deviation = 0.97
    16:13:52.644 INFO VariantDataManager - ReadPosRankSum: mean = 0.15 standard deviation = 1.02
    16:13:52.712 INFO VariantDataManager - Annotation order is: [DP, MQ, QD, FS, ReadPosRankSum]
    16:13:52.718 INFO VariantDataManager - Training with 13094 variants after standard deviation thresholding.
    16:13:52.720 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
    16:13:52.940 INFO VariantRecalibratorEngine - Finished iteration 0.
    16:13:53.076 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.07580
    16:13:53.165 INFO VariantRecalibratorEngine - Finished iteration 10. Current change in mixture coefficients = 0.06460
    16:13:53.253 INFO VariantRecalibratorEngine - Finished iteration 15. Current change in mixture coefficients = 0.01033
    16:13:53.341 INFO VariantRecalibratorEngine - Finished iteration 20. Current change in mixture coefficients = 0.01247
    16:13:53.414 INFO VariantRecalibratorEngine - Convergence after 24 iterations!
    16:13:53.434 INFO VariantRecalibratorEngine - Evaluating full set of 45108 variants...
    16:13:54.081 INFO VariantDataManager - Selected worst 1354 scoring variants --> variants with LOD <= -5.0000.
    16:13:54.081 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
    16:13:54.098 INFO VariantRecalibratorEngine - Finished iteration 0.
    16:13:54.104 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.02325
    16:13:54.110 INFO VariantRecalibratorEngine - Finished iteration 10. Current change in mixture coefficients = 0.00563
    16:13:54.117 INFO VariantRecalibratorEngine - Finished iteration 15. Current change in mixture coefficients = 0.00108
    16:13:54.117 INFO VariantRecalibratorEngine - Convergence after 15 iterations!
    16:13:54.127 INFO VariantRecalibratorEngine - Evaluating full set of 45108 variants...
    16:13:54.784 INFO TrancheManager - Finding 6 tranches for 45108 variants
    16:13:54.809 INFO TrancheManager - TruthSensitivityTranche threshold 60.00 => selection metric threshold 0.400
    16:13:54.821 INFO TrancheManager - Found tranche for 60.000: 0.400 threshold starting with variant 38815; running score is 0.400
    16:13:54.821 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=60.00 minVQSLod=19.2431 known=(6221 @ 2.0883) novel=(72 @ 0.5000) truthSites(6689 accessible, 4013 called), name=anonymous]
    16:13:54.821 INFO TrancheManager - TruthSensitivityTranche threshold 55.90 => selection metric threshold 0.441
    16:13:54.827 INFO TrancheManager - Found tranche for 55.900: 0.441 threshold starting with variant 39577; running score is 0.441
    16:13:54.827 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=55.90 minVQSLod=19.4768 known=(5517 @ 2.1017) novel=(14 @ 1.3333) truthSites(6689 accessible, 3739 called), name=anonymous]
    16:13:54.827 INFO TrancheManager - TruthSensitivityTranche threshold 55.80 => selection metric threshold 0.442
    16:13:54.830 INFO TrancheManager - Found tranche for 55.800: 0.442 threshold starting with variant 39592; running score is 0.442
    16:13:54.830 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=55.80 minVQSLod=19.4814 known=(5503 @ 2.1025) novel=(13 @ 1.6000) truthSites(6689 accessible, 3732 called), name=anonymous]
    16:13:54.830 INFO TrancheManager - TruthSensitivityTranche threshold 55.70 => selection metric threshold 0.443
    16:13:54.833 INFO TrancheManager - Found tranche for 55.700: 0.443 threshold starting with variant 39607; running score is 0.443
    16:13:54.833 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=55.70 minVQSLod=19.4869 known=(5489 @ 2.1034) novel=(12 @ 2.0000) truthSites(6689 accessible, 3725 called), name=anonymous]
    16:13:54.833 INFO TrancheManager - TruthSensitivityTranche threshold 55.50 => selection metric threshold 0.445
    16:13:54.836 INFO TrancheManager - Found tranche for 55.500: 0.445 threshold starting with variant 39643; running score is 0.445
    16:13:54.836 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=55.50 minVQSLod=19.4990 known=(5455 @ 2.1041) novel=(10 @ 2.3333) truthSites(6689 accessible, 3712 called), name=anonymous]
    16:13:54.836 INFO TrancheManager - TruthSensitivityTranche threshold 55.00 => selection metric threshold 0.450
    16:13:54.839 INFO TrancheManager - Found tranche for 55.000: 0.450 threshold starting with variant 39744; running score is 0.450
    16:13:54.839 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=55.00 minVQSLod=19.5381 known=(5359 @ 2.0886) novel=(5 @ 5.0000) truthSites(6689 accessible, 3678 called), name=anonymous]
    16:13:54.840 INFO VariantRecalibrator - Writing out recalibration table...
    16:13:55.390 INFO VariantRecalibrator - Writing out visualization Rscript file...
    16:13:55.406 INFO VariantRecalibrator - Building DP x MQ plot...
    16:13:55.410 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:55.856 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:56.079 INFO VariantRecalibrator - Building DP x QD plot...
    16:13:56.082 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:56.424 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:56.700 INFO VariantRecalibrator - Building DP x FS plot...
    16:13:56.703 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:57.044 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:57.255 INFO VariantRecalibrator - Building DP x ReadPosRankSum plot...
    16:13:57.257 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:57.608 INFO VariantRecalibratorEngine - Evaluating full set of 3660 variants...
    16:13:57.810 INFO VariantRecalibrator - Building MQ x QD plot...
    16:13:57.811 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:13:58.207 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:13:58.395 INFO VariantRecalibrator - Building MQ x FS plot...
    16:13:58.396 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:13:58.707 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:13:58.919 INFO VariantRecalibrator - Building MQ x ReadPosRankSum plot...
    16:13:58.921 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:13:59.846 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:00.017 INFO VariantRecalibrator - Building QD x FS plot...
    16:14:00.018 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:00.301 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:00.485 INFO VariantRecalibrator - Building QD x ReadPosRankSum plot...
    16:14:00.485 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:00.768 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:00.944 INFO VariantRecalibrator - Building FS x ReadPosRankSum plot...
    16:14:00.944 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:01.223 INFO VariantRecalibratorEngine - Evaluating full set of 3721 variants...
    16:14:01.405 INFO VariantRecalibrator - Executing: Rscript /mnt1/fastq/hs38dr/tmp_Merged_snps.plot.R
    16:14:12.480 INFO VariantRecalibrator - Executing: Rscript (resource)org/broadinstitute/hellbender/tools/walkers/vqsr/plot_Tranches.R /mnt1/fastq/hs38dr/tmp_sitesonly.snps.tranches 2.15
  • xiuczxiucz Member
    edited April 8

    @bhanuGandham
    The same problem I have met, it worked well with GATK4.0.11.0 and previous versions. But in the new version GATK4.1.0.0, it reported these warnings.

  • wlaiwlai Member

    @bhanuGandham i have also tried annotate them using parenthesis just so you know that it doesnt work for me too, wondering if @xiucz has tried with parentheses yet.
    "(QD < 2.0) || (FS > 60.0) || (MQ < 40.0) || (MQRankSum < -12.5) || ..."

  • xiuczxiucz Member

    @wlai sorry, I have not tried with that, but I think it is more meaningful and clear to split the regex with filter-name annotation labeled in the FILTER column.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @xiucz and @wlai

    The WARN statements simply let you know that the annotation is missing at the site. There is no need to to worry about them, as sometimes an annotation cannot be calculated for a site. Have a look at this thread: https://gatkforums.broadinstitute.org/gatk/discussion/2334/undefined-variable-variantfiltration

  • wlaiwlai Member

    @bhanuGandham
    Do you mean i should put less important annotation at the first position so that it is alright to skip it?

Sign In or Register to comment.