Service notice: Several of our team members are on vacation so service will be slow through at least July 13th, possibly longer depending on how much backlog accumulates during that time. This means that for a while it may take us more time than usual to answer your questions. Thank you for your patience.

Unexpected result for VQSR with non-model organism (Plasmodium falciparum - Malaria parasite).

Hi,

I am variant calling wgs (~best practices; group calling for g.vcf) in Plasmodium falciparum and I would like to filter the call data using VariantRecalibrator VQSR, as opposed to setting hard filters.

I have obtained some plasmodium training/truth/known score calibration data sets looking at crosses between different strains of Plasmodium - these are from the Sanger group (UK).

However, when I set up my data to be filtered with VQSR and run it, I find that it almost gives the opposite result to the one I am expecting.

My test (positive) SNPs fail the filter (aka PASS right through) and many highly questionable variants with low coverage are flagged as Significant by the filter. This is wrong.

Here is an example of a positive control snp failing:

Pf3D7_06_v3 131212 . C T 36619.8 PASS AC=14;AF=0.137;AN=102;BaseQRankSum=1.22;ClippingRankSum=0.306;DP=4115;ExcessHet=0.0000;FS=0.000;InbreedingCoeff=0.7990;MLEAC=14;MLEAF=0.137;MQ=60.00;MQRankSum=1.69;QD=29.34;ReadPosRankSum=1.34;SOR=0.565;VQSLOD=4.69;culprit=MQ GT:AD:DP:GQ:PL 0/0:46,0:46:99:0,111,1665 0/0:51,0:51:99:0,115,1794 0/0:54,0:54:99:0,116,1800 ./.:0,0:0 ./.:0,0:0 0/0:46,0:46:99:0,102,1539 1/1:0,58:58:99:2380,174,0 0/0:77,0:77:99:0,120,1800 ./. 0/0:58,0:58:99:0,112,1800 0/0:10,0:10:0:0,0,225 0/0:48,0:48:99:0,120,1800 0/0:45,0:45:99:0,119,1525 0/0:67,0:67:99:0,120,1800 ./.:0,0:0 0/0:26,0:26:72:0,72,890 0/0:135,0:135:99:0,120,1800 0/0:48,0:48:99:0,108,1362 ./.:0,0:0 0/0:110,0:110:99:0,120,1800 0/0:66,0:66:99:0,120,1800 0/0:57,0:57:99:0,120,1800 0/0:59,0:59:99:0,103,1800 0/0:70,0:70:99:0,118,1800 0/0:86,0:86:99:0,120,1800 0/0:67,0:67:99:0,118,1800 0/0:60,0:60:99:0,120,1800 0/0:54,0:54:99:0,114,1800 0/0:109,0:109:99:0,120,1800 0/0:65,0:65:99:0,105,1800 0/0:84,0:84:99:0,120,1800 0/0:66,0:66:99:0,101,1800 0/0:60,0:60:99:0,120,1800 0/0:66,0:66:99:0,103,1800 0/0:105,0:105:99:0,120,1800 0/0:60,0:60:99:0,113,1800 0/0:65,0:65:99:0,111,1800 0/0:50,0:50:99:0,117,1755 0/0:60,0:60:99:0,120,1800 0/1:109,59:168:99:1551,0,3432 0/0:114,0:114:99:0,120,1800 1/1:0,160:160:99:5769,480,0 1/1:0,150:150:99:5256,451,0 0/1:72,101:173:99:3043,0,1914 0/0:86,0:86:99:0,120,1800 1/1:0,180:180:99:6180,541,0 1/1:0,166:166:99:5470,495,0 1/1:0,193:193:99:7154,580,0 0/0:86,0:86:99:0,120,1800 0/0:62,0:62:99:0,120,1800 0/0:48,0:48:99:0,120,1800 0/0:59,0:59:99:0,109,1800 0/0:93,0:93:99:0,120,1800 0/0:51,0:51:99:0,118,1800 0/0:48,0:48:99:0,120,1800 0/0:61,0:61:99:0,106,1800

Here is a background call (false positive) being picked up as significant:

Pf3D7_06_v3 130463 . G T 120.41 VQSRTrancheSNP99.90to100.00 AC=2;AF=0.020;AN=102;BaseQRankSum=0.317;ClippingRankSum=0.917;DP=3763;ExcessHet=3.1387;FS=4.969;InbreedingCoeff=-0.0251;MLEAC=2;MLEAF=0.020;MQ=60.00;MQRankSum=1.35;NEGATIVE_TRAIN_SITE;QD=0.34;ReadPosRankSum=1.11;SOR=0.796;VQSLOD=-5.004e+00;culprit=QD GT:AD:DP:GQ:PL 0/0:58,0:58:99:0,118,1800 0/0:60,0:60:99:0,111,1636 0/0:64,0:64:99:0,120,1800 ./.:0,0:0 ./.:0,0:0 0/0:57,0:57:99:0,100,1431 0/0:46,0:46:99:0,105,1575 0/0:72,0:72:99:0,120,1800 ./. 0/0:48,0:48:99:0,111,1665 0/0:17,0:17:0:0,0,385 0/0:68,0:68:99:0,113,1800 0/0:73,0:73:99:0,120,1800 0/0:48,0:48:99:0,120,1800 ./.:0,0:0 0/0:21,0:21:60:0,60,900 0/0:183,0:183:99:0,120,1800 0/0:47,0:47:99:0,105,1436 ./.:0,0:0 0/0:40,0:40:99:0,100,1114 0/0:49,0:49:99:0,119,1800 0/0:62,0:62:99:0,102,1800 0/0:83,0:83:99:0,120,1800 0/0:93,0:93:99:0,119,1800 0/0:60,0:60:99:0,104,1800 0/0:45,0:45:99:0,114,1065 0/0:107,0:107:99:0,120,1800 0/0:76,0:76:99:0,116,1800 0/0:68,0:68:99:0,114,1800 0/0:73,0:73:99:0,102,1800 0/0:67,0:67:99:0,116,1800 0/1:137,17:154:81:81,0,3526 0/0:75,0:75:99:0,99,1485 0/0:76,0:76:99:0,113,1800 0/1:182,23:205:93:93,0,5061 0/0:64,0:64:99:0,120,1800 0/0:65,0:65:99:0,111,1800 0/0:85,0:85:99:0,115,1800 0/0:82,0:82:99:0,101,1800 0/0:70,0:70:99:0,120,1800 0/0:84,0:84:99:0,120,1800 0/0:66,0:66:99:0,120,1800 0/0:52,0:52:99:0,116,1493 0/0:75,0:75:99:0,120,1800 0/0:127,0:127:99:0,120,1800 0/0:68,0:68:99:0,120,1800 0/0:70,0:70:99:0,120,1800 0/0:55,0:55:99:0,119,1800 0/0:61,0:61:99:0,120,1800 0/0:73,0:73:99:0,119,1800 0/0:78,0:78:99:0,120,1800 0/0:59,0:59:99:0,109,1800 0/0:77,0:77:99:0,120,1800 0/0:102,0:102:99:0,120,1800 0/0:75,0:75:99:0,120,1800 0/0:74,0:74:99:0,120,1800

I would be very grateful if you have any suggestions as to how I should set this up - -or trouble shoot what I already have.

I could make another truth-calibration set since I have >300 snps and many indels which have been sequence verified. Is there a guide on constructing a vqsr calibration resource?

Many thanks for any information,
Regards,
Roy Williams

This is the command I used in gatk4.0

REFSNP1="$KS/3d7_hb3.combined.final.vcf.gz"
REFSNP2="$KS/7g8_gb4.combined.final.vcf.gz"
REFSNP3="$KS/hb3_dd2.combined.final.vcf.gz"

[[email protected] pf_variant_calling_pipeline-master]$ gatk VariantRecalibrator \

-R $REF \
-V $input \
--resource cross1,training=true,truth=true,known=false:$REFSNP1 \
--resource cross2,training=true,truth=true,known=false:$REFSNP2 \
--resource cross3,training=true,truth=true,known=false:$REFSNP3 \
--max-gaussians 4 \
-an QD \
-an MQ \
-an FS \
-an SOR \
-an DP \
-mode SNP \
-mq-cap 70 \
--output $output.snp.recal \
--tranches-file $output.snp.tranches

Using GATK jar /opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar
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=1 -jar /opt/applicationsiantRecalibrator -R /storage/NFS/GENOME_RESOURCES/pf/plasmodb-download/Pfalciparum3D7/fasta/data/PlasmoDB-30_Pfalciparum3D7_Genome.fasta -V /storage/extra-home/rmw002/PROJECTS_INFORMATICect/plasmodium_snp_analysis/output/PF_6sets_no_p_viv/combine_GVCFs/OUTFILE_1.snps.indels.g.vcf --resource cross1,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_ant_calling_pipeline-master/data/known_sites//3d7_hb3.combined.final.vcf.gz --resource cross2,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODe-master/data/known_sites//7g8_gb4.combined.final.vcf.gz --resource cross3,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2_sites//hb3_dd2.combined.final.vcf.gz --max-gaussians 4 -an QD -an MQ -an FS -an SOR -an DP -mode SNP -mq-cap 70 --output OUTFILE_1_VQSR.snp.recal --tranches-file OUTFILE_1_VQSR.snp.tran
17:08:43.858 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar!/com/intel/gkl/native/libgkl_compressio
17:08:43.952 INFO VariantRecalibrator - ------------------------------------------------------------
17:08:43.953 INFO VariantRecalibrator - The Genome Analysis Toolkit (GATK) v4.0.1.2
17:08:43.953 INFO VariantRecalibrator - For support and documentation go to https://software.broadinstitute.org/gatk/
17:08:43.953 INFO VariantRecalibrator - Executing as [email protected] on Linux v3.10.0-327.36.3.el7.x86_64 amd64
17:08:43.953 INFO VariantRecalibrator - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_65-b17
17:08:43.953 INFO VariantRecalibrator - Start Date/Time: February 21, 2018 5:08:43 PM PST
17:08:43.953 INFO VariantRecalibrator - ------------------------------------------------------------
17:08:43.953 INFO VariantRecalibrator - ------------------------------------------------------------
17:08:43.954 INFO VariantRecalibrator - HTSJDK Version: 2.14.1
17:08:43.954 INFO VariantRecalibrator - Picard Version: 2.17.2
17:08:43.954 INFO VariantRecalibrator - HTSJDK Defaults.COMPRESSION_LEVEL : 1
17:08:43.954 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:08:43.954 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:08:43.954 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:08:43.954 INFO VariantRecalibrator - Deflater: IntelDeflater
17:08:43.954 INFO VariantRecalibrator - Inflater: IntelInflater
17:08:43.954 INFO VariantRecalibrator - GCS max retries/reopens: 20
17:08:43.954 INFO VariantRecalibrator - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
17:08:43.954 INFO VariantRecalibrator - Initializing engine
17:08:44.596 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:08:44.617 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:08:44.627 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:08:44.635 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/root_rstudio_project/plasmodium_CFs/OUTFILE_1.snps.indels.g.vcf
17:08:44.661 INFO VariantRecalibrator - Done initializing engine
17:08:44.663 INFO TrainingSet - Found cross1 track: Known = false Training = true Truth = true Prior = Q0.0
17:08:44.663 INFO TrainingSet - Found cross2 track: Known = false Training = true Truth = true Prior = Q0.0
17:08:44.663 INFO TrainingSet - Found cross3 track: Known = false Training = true Truth = true Prior = Q0.0
17:08:44.667 WARN GATKVariantContextUtils - Can't determine output variant file format from output file extension "recal". Defaulting to VCF.
17:08:44.683 INFO ProgressMeter - Starting traversal
17:08:44.683 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:08:54.046 INFO ProgressMeter - Pf3D7_14_v3:3283328 0.2 263809 1690722.1
17:08:54.046 INFO ProgressMeter - Traversal complete. Processed 263809 total variants in 0.2 minutes.
17:08:54.062 INFO VariantDataManager - QD: mean = 29.99 standard deviation = 4.26
17:08:54.081 INFO VariantDataManager - MQ: mean = 1.56 standard deviation = 0.57
17:08:54.098 INFO VariantDataManager - FS: mean = 1.25 standard deviation = 5.54
17:08:54.107 INFO VariantDataManager - SOR: mean = 0.96 standard deviation = 0.58
17:08:54.116 INFO VariantDataManager - DP: mean = 2521.37 standard deviation = 1188.83
17:08:54.378 INFO VariantDataManager - Annotations are now ordered by their information content: [DP, QD, MQ, SOR, FS]
17:08:54.387 INFO VariantDataManager - Training with 13738 variants after standard deviation thresholding.
17:08:54.390 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
17:08:54.753 INFO VariantRecalibratorEngine - Finished iteration 0.
17:08:54.896 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.19610
17:08:54.996 INFO VariantRecalibratorEngine - Finished iteration 10. Current change in mixture coefficients = 0.03783
17:08:55.097 INFO VariantRecalibratorEngine - Finished iteration 15. Current change in mixture coefficients = 0.03810
17:08:55.199 INFO VariantRecalibratorEngine - Finished iteration 20. Current change in mixture coefficients = 0.01644
17:08:55.301 INFO VariantRecalibratorEngine - Finished iteration 25. Current change in mixture coefficients = 0.00708
17:08:55.401 INFO VariantRecalibratorEngine - Finished iteration 30. Current change in mixture coefficients = 0.00391
17:08:55.496 INFO VariantRecalibratorEngine - Finished iteration 35. Current change in mixture coefficients = 0.00230
17:08:55.534 INFO VariantRecalibratorEngine - Convergence after 37 iterations!
17:08:55.560 INFO VariantRecalibratorEngine - Evaluating full set of 148347 variants...
17:08:55.782 INFO VariantDataManager - Selected worst 61509 scoring variants --> variants with LOD <= -5.0000.
17:08:55.782 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
17:08:57.140 INFO VariantRecalibratorEngine - Finished iteration 0.
17:08:57.519 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.03811
17:08:57.897 INFO VariantRecalibratorEngine - Finished iteration 10. Current change in mixture coefficients = 0.02126
17:08:58.285 INFO VariantRecalibratorEngine - Finished iteration 15. Current change in mixture coefficients = 0.01184
17:08:58.670 INFO VariantRecalibratorEngine - Finished iteration 20. Current change in mixture coefficients = 0.00539
17:08:59.049 INFO VariantRecalibratorEngine - Finished iteration 25. Current change in mixture coefficients = 0.00180
17:08:59.049 INFO VariantRecalibratorEngine - Convergence after 25 iterations!
17:08:59.128 INFO VariantRecalibratorEngine - Evaluating full set of 148347 variants...
17:08:59.912 INFO TrancheManager - Finding 4 tranches for 148347 variants
17:09:00.003 INFO TrancheManager - TruthSensitivityTranche threshold 100.00 => selection metric threshold 0.000
17:09:00.044 INFO TrancheManager - Found tranche for 100.000: 0.000 threshold starting with variant 0; running score is 0.000
17:09:00.044 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=100.00 minVQSLod=-1317.6093 known=(0 @ 0.0000) novel=(148347 @ 0.7033) truonymous]
17:09:00.045 INFO TrancheManager - TruthSensitivityTranche threshold 99.90 => selection metric threshold 0.001
17:09:00.079 INFO TrancheManager - Found tranche for 99.900: 0.001 threshold starting with variant 8045; running score is 0.001
17:09:00.079 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=99.90 minVQSLod=-11.9390 known=(0 @ 0.0000) novel=(140302 @ 0.7067) truthSmous]
17:09:00.082 INFO TrancheManager - TruthSensitivityTranche threshold 99.00 => selection metric threshold 0.010
17:09:00.098 INFO TrancheManager - Found tranche for 99.000: 0.010 threshold starting with variant 49975; running score is 0.010
17:09:00.099 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2376 known=(0 @ 0.0000) novel=(98372 @ 0.8927) truthSitus]
17:09:00.099 INFO TrancheManager - TruthSensitivityTranche threshold 90.00 => selection metric threshold 0.100
17:09:00.112 INFO TrancheManager - Found tranche for 90.000: 0.100 threshold starting with variant 99056; running score is 0.100
17:09:00.112 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=90.00 minVQSLod=1.6878 known=(0 @ 0.0000) novel=(49291 @ 0.9689) truthSites]
17:09:00.113 INFO VariantRecalibrator - Writing out recalibration table...
17:09:01.162 INFO VariantRecalibrator - Shutting down engine
[February 21, 2018 5:09:01 PM PST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator done. Elapsed time: 0.29 minutes.
Runtime.totalMemory()=3971481600
Tool returned:
true
[[email protected] pf_variant_calling_pipeline-master]$ #
[[email protected] pf_variant_calling_pipeline-master]$ gatk VariantRecalibrator \

-R $REF \
-V $input \
--resource cross1,training=true,truth=true,known=false:$REFSNP1 \
--resource cross2,training=true,truth=true,known=false:$REFSNP2 \
--resource cross3,training=true,truth=true,known=false:$REFSNP3 \
--max-gaussians 4 \
-an QD \
-an MQ \
-an FS \
-an SOR \
-an DP \
-mode INDEL \
-mq-cap 70 \
--output $output.indel.recal \
--tranches-file $output.indel.tranches

Using GATK jar /opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar
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=1 -jar /opt/applicationsiantRecalibrator -R /storage/NFS/GENOME_RESOURCES/pf/plasmodb-download/Pfalciparum3D7/fasta/data/PlasmoDB-30_Pfalciparum3D7_Genome.fasta -V /storage/extra-home/rmw002/PROJECTS_INFORMATICect/plasmodium_snp_analysis/output/PF_6sets_no_p_viv/combine_GVCFs/OUTFILE_1.snps.indels.g.vcf --resource cross1,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_ant_calling_pipeline-master/data/known_sites//3d7_hb3.combined.final.vcf.gz --resource cross2,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODe-master/data/known_sites//7g8_gb4.combined.final.vcf.gz --resource cross3,training=true,truth=true,known=false:/storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2_sites//hb3_dd2.combined.final.vcf.gz --max-gaussians 4 -an QD -an MQ -an FS -an SOR -an DP -mode INDEL -mq-cap 70 --output OUTFILE_1_VQSR.indel.recal --tranches-file OUTFILE_1_VQSR.inde
17:16:42.523 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar!/com/intel/gkl/native/libgkl_compressio
17:16:42.619 INFO VariantRecalibrator - ------------------------------------------------------------
17:16:42.620 INFO VariantRecalibrator - The Genome Analysis Toolkit (GATK) v4.0.1.2
17:16:42.620 INFO VariantRecalibrator - For support and documentation go to https://software.broadinstitute.org/gatk/
17:16:42.620 INFO VariantRecalibrator - Executing as [email protected] on Linux v3.10.0-327.36.3.el7.x86_64 amd64
17:16:42.620 INFO VariantRecalibrator - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_65-b17
17:16:42.620 INFO VariantRecalibrator - Start Date/Time: February 21, 2018 5:16:42 PM PST
17:16:42.620 INFO VariantRecalibrator - ------------------------------------------------------------
17:16:42.620 INFO VariantRecalibrator - ------------------------------------------------------------
17:16:42.621 INFO VariantRecalibrator - HTSJDK Version: 2.14.1
17:16:42.621 INFO VariantRecalibrator - Picard Version: 2.17.2
17:16:42.621 INFO VariantRecalibrator - HTSJDK Defaults.COMPRESSION_LEVEL : 1
17:16:42.621 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:16:42.621 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:16:42.621 INFO VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:16:42.621 INFO VariantRecalibrator - Deflater: IntelDeflater
17:16:42.621 INFO VariantRecalibrator - Inflater: IntelInflater
17:16:42.621 INFO VariantRecalibrator - GCS max retries/reopens: 20
17:16:42.621 INFO VariantRecalibrator - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
17:16:42.621 INFO VariantRecalibrator - Initializing engine
17:16:43.034 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:16:43.055 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:16:43.065 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mastz
17:16:43.074 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/root_rstudio_project/plasmodium_CFs/OUTFILE_1.snps.indels.g.vcf
17:16:43.101 INFO VariantRecalibrator - Done initializing engine
17:16:43.102 INFO TrainingSet - Found cross1 track: Known = false Training = true Truth = true Prior = Q0.0
17:16:43.102 INFO TrainingSet - Found cross2 track: Known = false Training = true Truth = true Prior = Q0.0
17:16:43.103 INFO TrainingSet - Found cross3 track: Known = false Training = true Truth = true Prior = Q0.0
17:16:43.106 WARN GATKVariantContextUtils - Can't determine output variant file format from output file extension "recal". Defaulting to VCF.
17:16:43.122 INFO ProgressMeter - Starting traversal
17:16:43.122 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:16:52.707 INFO ProgressMeter - Pf3D7_14_v3:3283328 0.2 263809 1651558.8
17:16:52.707 INFO ProgressMeter - Traversal complete. Processed 263809 total variants in 0.2 minutes.
17:16:52.726 INFO VariantDataManager - QD: mean = 24.83 standard deviation = 8.97
17:16:52.742 INFO VariantDataManager - MQ: mean = 1.55 standard deviation = 0.51
17:16:52.756 INFO VariantDataManager - FS: mean = 3.17 standard deviation = 6.82
17:16:52.770 INFO VariantDataManager - SOR: mean = 0.94 standard deviation = 0.63
17:16:52.781 INFO VariantDataManager - DP: mean = 1511.44 standard deviation = 1023.70
17:16:52.994 INFO VariantDataManager - Annotations are now ordered by their information content: [DP, QD, FS, MQ, SOR]
17:16:53.004 INFO VariantDataManager - Training with 30537 variants after standard deviation thresholding.
17:16:53.007 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
17:16:53.583 INFO VariantRecalibratorEngine - Finished iteration 0.
17:16:53.919 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.21947
17:16:54.133 INFO VariantRecalibratorEngine - Finished iteration 10. Current change in mixture coefficients = 0.09436
17:16:54.349 INFO VariantRecalibratorEngine - Finished iteration 15. Current change in mixture coefficients = 0.05508
17:16:54.569 INFO VariantRecalibratorEngine - Finished iteration 20. Current change in mixture coefficients = 0.03812
17:16:54.788 INFO VariantRecalibratorEngine - Finished iteration 25. Current change in mixture coefficients = 0.03574
17:16:55.021 INFO VariantRecalibratorEngine - Finished iteration 30. Current change in mixture coefficients = 0.09227
17:16:55.267 INFO VariantRecalibratorEngine - Finished iteration 35. Current change in mixture coefficients = 0.01469
17:16:55.509 INFO VariantRecalibratorEngine - Finished iteration 40. Current change in mixture coefficients = 0.00464
17:16:55.750 INFO VariantRecalibratorEngine - Finished iteration 45. Current change in mixture coefficients = 0.00220
17:16:55.799 INFO VariantRecalibratorEngine - Convergence after 46 iterations!
17:16:55.849 INFO VariantRecalibratorEngine - Evaluating full set of 115462 variants...
17:16:55.996 INFO VariantDataManager - Selected worst 16780 scoring variants --> variants with LOD <= -5.0000.
17:16:55.996 INFO GaussianMixtureModel - Initializing model with 100 k-means iterations...
17:16:56.207 INFO VariantRecalibratorEngine - Finished iteration 0.
17:16:56.260 INFO VariantRecalibratorEngine - Finished iteration 5. Current change in mixture coefficients = 0.00104
17:16:56.260 INFO VariantRecalibratorEngine - Convergence after 5 iterations!
17:16:56.271 INFO VariantRecalibratorEngine - Evaluating full set of 115462 variants...
17:16:56.794 INFO TrancheManager - Finding 4 tranches for 115462 variants
17:16:56.862 INFO TrancheManager - TruthSensitivityTranche threshold 100.00 => selection metric threshold 0.000
17:16:56.891 INFO TrancheManager - Found tranche for 100.000: 0.000 threshold starting with variant 0; running score is 0.000
17:16:56.891 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=100.00 minVQSLod=-479.2352 known=(0 @ 0.0000) novel=(115462 @ 0.0000) trutnymous]
17:16:56.891 INFO TrancheManager - TruthSensitivityTranche threshold 99.90 => selection metric threshold 0.001
17:16:56.908 INFO TrancheManager - Found tranche for 99.900: 0.001 threshold starting with variant 1556; running score is 0.001
17:16:56.908 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=99.90 minVQSLod=-8.3698 known=(0 @ 0.0000) novel=(113906 @ 0.0000) truthSious]
17:16:56.908 INFO TrancheManager - TruthSensitivityTranche threshold 99.00 => selection metric threshold 0.010
17:16:56.919 INFO TrancheManager - Found tranche for 99.000: 0.010 threshold starting with variant 11164; running score is 0.010
17:16:56.922 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2627 known=(0 @ 0.0000) novel=(104298 @ 0.0000) truthSious]
17:16:56.922 INFO TrancheManager - TruthSensitivityTranche threshold 90.00 => selection metric threshold 0.100
17:16:56.933 INFO TrancheManager - Found tranche for 90.000: 0.100 threshold starting with variant 37427; running score is 0.100
17:16:56.933 INFO TrancheManager - TruthSensitivityTranche is TruthSensitivityTranche targetTruthSensitivity=90.00 minVQSLod=0.6190 known=(0 @ 0.0000) novel=(78035 @ 0.0000) truthSites]
17:16:56.934 INFO VariantRecalibrator - Writing out recalibration table...
17:16:57.840 INFO VariantRecalibrator - Tranches plot will not be generated since we are running in INDEL mode
17:16:57.853 INFO VariantRecalibrator - Shutting down engine
[February 21, 2018 5:16:57 PM PST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator done. Elapsed time: 0.26 minutes.
Runtime.totalMemory()=4442816512
Tool returned:
true
[[email protected] pf_variant_calling_pipeline-master]$ #SNP
[[email protected] pf_variant_calling_pipeline-master]$ gatk ApplyVQSR \

-R $REF \
-V $input \
-O $output.snp.vcf \
-ts-filter-level 99.0 \
--tranches-file $output.snp.tranches \
--recal-file $output.snp.recal \
-mode SNP

Using GATK jar /opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar
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=1 -jar /opt/applicationslyVQSR -R /storage/NFS/GENOME_RESOURCES/pf/plasmodb-download/Pfalciparum3D7/fasta/data/PlasmoDB-30_Pfalciparum3D7_Genome.fasta -V /storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIdium_snp_analysis/output/PF_6sets_no_p_viv/combine_GVCFs/OUTFILE_1.snps.indels.g.vcf -O OUTFILE_1_VQSR.snp.vcf -ts-filter-level 99.0 --tranches-file OUTFILE_1_VQSR.snp.tranches --recal-f
17:17:54.203 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar!/com/intel/gkl/native/libgkl_compressio
17:17:54.299 INFO ApplyVQSR - ------------------------------------------------------------
17:17:54.299 INFO ApplyVQSR - The Genome Analysis Toolkit (GATK) v4.0.1.2
17:17:54.300 INFO ApplyVQSR - For support and documentation go to https://software.broadinstitute.org/gatk/
17:17:54.300 INFO ApplyVQSR - Executing as [email protected] on Linux v3.10.0-327.36.3.el7.x86_64 amd64
17:17:54.300 INFO ApplyVQSR - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_65-b17
17:17:54.300 INFO ApplyVQSR - Start Date/Time: February 21, 2018 5:17:54 PM PST
17:17:54.300 INFO ApplyVQSR - ------------------------------------------------------------
17:17:54.300 INFO ApplyVQSR - ------------------------------------------------------------
17:17:54.301 INFO ApplyVQSR - HTSJDK Version: 2.14.1
17:17:54.301 INFO ApplyVQSR - Picard Version: 2.17.2
17:17:54.301 INFO ApplyVQSR - HTSJDK Defaults.COMPRESSION_LEVEL : 1
17:17:54.301 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:17:54.301 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:17:54.301 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:17:54.301 INFO ApplyVQSR - Deflater: IntelDeflater
17:17:54.301 INFO ApplyVQSR - Inflater: IntelInflater
17:17:54.301 INFO ApplyVQSR - GCS max retries/reopens: 20
17:17:54.301 INFO ApplyVQSR - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
17:17:54.301 INFO ApplyVQSR - Initializing engine
17:17:54.942 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mast
17:17:54.960 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/root_rstudio_project/plasmodium_CFs/OUTFILE_1.snps.indels.g.vcf
17:17:54.986 INFO ApplyVQSR - Done initializing engine
17:17:54.988 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=90.00 minVQSLod=1.6878 known=(0 @ 0.0000) novel=(49291 @ 0.9689) truthSites(13760 accessible, 1
17:17:54.988 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2376 known=(0 @ 0.0000) novel=(98372 @ 0.8927) truthSites(13760 accessible,
17:17:54.989 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=99.90 minVQSLod=-11.9390 known=(0 @ 0.0000) novel=(140302 @ 0.7067) truthSites(13760 accessible]
17:17:54.989 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=100.00 minVQSLod=-1317.6093 known=(0 @ 0.0000) novel=(148347 @ 0.7033) truthSites(13760 accessi0.00]
17:17:55.004 INFO ApplyVQSR - Keeping all variants in tranche TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2376 known=(0 @ 0.0000) novel=(98372 @ 0.8927) truthSitesheSNP90.00to99.00]
17:17:55.087 INFO ProgressMeter - Starting traversal
17:17:55.087 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:18:03.084 INFO ProgressMeter - Pf3D7_14_v3:3283328 0.1 263809 1979557.3
17:18:03.084 INFO ProgressMeter - Traversal complete. Processed 263809 total variants in 0.1 minutes.
17:18:03.227 INFO ApplyVQSR - Shutting down engine
[February 21, 2018 5:18:03 PM PST] org.broadinstitute.hellbender.tools.walkers.vqsr.ApplyVQSR done. Elapsed time: 0.15 minutes.
Runtime.totalMemory()=3329228800
[[email protected] pf_variant_calling_pipeline-master]$ #
[[email protected] pf_variant_calling_pipeline-master]$ #INDEL
[[email protected] pf_variant_calling_pipeline-master]$ gatk ApplyVQSR \

-R $REF \
-V $input \
-O $output.indels.vcf \
--ts-filter-level 99.0 \
--tranches-file $output.indel.tranches \
--recal-file $output.indel.recal \
-mode INDEL

Using GATK jar /opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar
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=1 -jar /opt/applicationslyVQSR -R /storage/NFS/GENOME_RESOURCES/pf/plasmodb-download/Pfalciparum3D7/fasta/data/PlasmoDB-30_Pfalciparum3D7_Genome.fasta -V /storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIdium_snp_analysis/output/PF_6sets_no_p_viv/combine_GVCFs/OUTFILE_1.snps.indels.g.vcf -O OUTFILE_1_VQSR.indels.vcf --ts-filter-level 99.0 --tranches-file OUTFILE_1_VQSR.indel.tranches --r
17:18:06.051 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar!/com/intel/gkl/native/libgkl_compressio
17:18:06.144 INFO ApplyVQSR - ------------------------------------------------------------
17:18:06.144 INFO ApplyVQSR - The Genome Analysis Toolkit (GATK) v4.0.1.2
17:18:06.144 INFO ApplyVQSR - For support and documentation go to https://software.broadinstitute.org/gatk/
17:18:06.144 INFO ApplyVQSR - Executing as [email protected] on Linux v3.10.0-327.36.3.el7.x86_64 amd64
17:18:06.144 INFO ApplyVQSR - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_65-b17
17:18:06.145 INFO ApplyVQSR - Start Date/Time: February 21, 2018 5:18:06 PM PST
17:18:06.145 INFO ApplyVQSR - ------------------------------------------------------------
17:18:06.145 INFO ApplyVQSR - ------------------------------------------------------------
17:18:06.145 INFO ApplyVQSR - HTSJDK Version: 2.14.1
17:18:06.145 INFO ApplyVQSR - Picard Version: 2.17.2
17:18:06.145 INFO ApplyVQSR - HTSJDK Defaults.COMPRESSION_LEVEL : 1
17:18:06.145 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:18:06.145 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:18:06.145 INFO ApplyVQSR - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:18:06.146 INFO ApplyVQSR - Deflater: IntelDeflater
17:18:06.146 INFO ApplyVQSR - Inflater: IntelInflater
17:18:06.146 INFO ApplyVQSR - GCS max retries/reopens: 20
17:18:06.146 INFO ApplyVQSR - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
17:18:06.146 INFO ApplyVQSR - Initializing engine
17:18:06.806 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mast
17:18:06.833 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/root_rstudio_project/plasmodium_CFs/OUTFILE_1.snps.indels.g.vcf
17:18:06.864 INFO ApplyVQSR - Done initializing engine
17:18:06.867 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=90.00 minVQSLod=0.6190 known=(0 @ 0.0000) novel=(78035 @ 0.0000) truthSites(30572 accessible, 2
17:18:06.867 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2627 known=(0 @ 0.0000) novel=(104298 @ 0.0000) truthSites(30572 accessible,0]
17:18:06.867 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=99.90 minVQSLod=-8.3698 known=(0 @ 0.0000) novel=(113906 @ 0.0000) truthSites(30572 accessible,0]
17:18:06.868 INFO ApplyVQSR - Read tranche TruthSensitivityTranche targetTruthSensitivity=100.00 minVQSLod=-479.2352 known=(0 @ 0.0000) novel=(115462 @ 0.0000) truthSites(30572 accessib00.00]
17:18:06.884 INFO ApplyVQSR - Keeping all variants in tranche TruthSensitivityTranche targetTruthSensitivity=99.00 minVQSLod=-2.2627 known=(0 @ 0.0000) novel=(104298 @ 0.0000) truthSitecheINDEL90.00to99.00]
17:18:06.969 INFO ProgressMeter - Starting traversal
17:18:06.969 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:18:15.685 INFO ProgressMeter - Pf3D7_14_v3:3283328 0.1 263809 1816241.0
17:18:15.685 INFO ProgressMeter - Traversal complete. Processed 263809 total variants in 0.1 minutes.
17:18:15.835 INFO ApplyVQSR - Shutting down engine
[February 21, 2018 5:18:15 PM PST] org.broadinstitute.hellbender.tools.walkers.vqsr.ApplyVQSR done. Elapsed time: 0.16 minutes.
Runtime.totalMemory()=4833411072
[[email protected] pf_variant_calling_pipeline-master]$ #
[[email protected] pf_variant_calling_pipeline-master]$ #SNP
[[email protected] pf_variant_calling_pipeline-master]$ gatk SelectVariants \

 -R $REF \
 -V OUTFILE_1_VQSR.snp.vcf \
 --select-type-to-include SNP \
 -O $output.snp.99.vqsr.vcf

Using GATK jar /opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar
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=1 -jar /opt/applicationsectVariants -R /storage/NFS/GENOME_RESOURCES/pf/plasmodb-download/Pfalciparum3D7/fasta/data/PlasmoDB-30_Pfalciparum3D7_Genome.fasta -V OUTFILE_1_VQSR.snp.vcf --select-type-to-include SNP
17:18:18.713 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/applications/gatk/4.0.1.2/gatk-package-4.0.1.2-local.jar!/com/intel/gkl/native/libgkl_compressio
17:18:18.807 INFO SelectVariants - ------------------------------------------------------------
17:18:18.808 INFO SelectVariants - The Genome Analysis Toolkit (GATK) v4.0.1.2
17:18:18.808 INFO SelectVariants - For support and documentation go to https://software.broadinstitute.org/gatk/
17:18:18.808 INFO SelectVariants - Executing as [email protected] on Linux v3.10.0-327.36.3.el7.x86_64 amd64
17:18:18.808 INFO SelectVariants - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_65-b17
17:18:18.808 INFO SelectVariants - Start Date/Time: February 21, 2018 5:18:18 PM PST
17:18:18.808 INFO SelectVariants - ------------------------------------------------------------
17:18:18.808 INFO SelectVariants - ------------------------------------------------------------
17:18:18.809 INFO SelectVariants - HTSJDK Version: 2.14.1
17:18:18.809 INFO SelectVariants - Picard Version: 2.17.2
17:18:18.809 INFO SelectVariants - HTSJDK Defaults.COMPRESSION_LEVEL : 1
17:18:18.809 INFO SelectVariants - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:18:18.809 INFO SelectVariants - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:18:18.809 INFO SelectVariants - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:18:18.809 INFO SelectVariants - Deflater: IntelDeflater
17:18:18.809 INFO SelectVariants - Inflater: IntelInflater
17:18:18.809 INFO SelectVariants - GCS max retries/reopens: 20
17:18:18.809 INFO SelectVariants - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
17:18:18.809 INFO SelectVariants - Initializing engine
17:18:19.307 INFO FeatureManager - Using codec VCFCodec to read file file:///storage/extra-home/rmw002/PROJECTS_INFORMATICS/PLASMODIUM/SNP_PIPELINE_v2.0/pf_variant_calling_pipeline-mast
17:18:19.336 INFO SelectVariants - Done initializing engine
17:18:19.410 INFO ProgressMeter - Starting traversal
17:18:19.411 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:18:29.422 INFO ProgressMeter - Pf3D7_08_v3:173442 0.2 66000 395604.4
17:18:39.033 INFO ProgressMeter - Pf3D7_14_v3:3291520 0.3 148347 453637.4
17:18:39.033 INFO ProgressMeter - Traversal complete. Processed 148347 total variants in 0.3 minutes.
17:18:39.093 INFO SelectVariants - Shutting down engine
[February 21, 2018 5:18:39 PM PST] org.broadinstitute.hellbender.tools.walkers.variantutils.SelectVariants done. Elapsed time: 0.34 minutes.
Runtime.totalMemory()=4976541696

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator
    edited February 28

    @royw
    Hi,

    I think you are confusing PASS and the VQSRTranche filter. PASS variants are the ones you keep/are called true positive. The ones that fail/are filtered out are labeled by the tranche they fall in. For example, in your second example record which is a false positive, notice that filter field contains VQSRTrancheSNP99.90to100.00. This means the VQSLOD score is in the bottom .1% compared to the VQSLOD scores of the sites in your truth set.

    I realize the 99.90 to 100.00 is confusing, but perhaps the VQSR presentation and hands on tutorial in the Presentations section will help.

    -Sheila

Sign In or Register to comment.