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.

Different results from identical GenotypeGVCFs runs when multithreading

JverlouwJverlouw Erasmus MC, RotterdamMember
edited April 2016 in Ask the GATK team

Dear GATK Devs,

While working on a new version of our human exome pipeline I ran into the problem that the output VCF files were not the same, even though the commands and samples run through it were exactly the same. Some of the quality scores like MQranksum or QD are slightly different, but this influences the end result quite a bit. Initially I thought this might be our older version of GATK (3.3), but I noticed the same effect when using 3.5.

I couldn't find anything about this on the forum or in the documentation. Is this a known bug by chance?

Below is the result of a diff command on the two resulting VCF files. Paths and sample IDs are censored, but are identical for both.

14c14
##GATKCommandLine.GenotypeGVCFs=<ID=GenotypeGVCFs,Version=3.5-0-g36282e4,Date="Thu Apr 21 11:55:39 BST 2016",Epoch=1461236139620,CommandLineOptions="analysis_type=GenotypeGVCFs input_file=[] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=[/Reference/SeqCap_44Mb_Exome.bed] excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/Reference/ucsc.hg19.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 static_quantized_quals=null round_down_quantized=false disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=true no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=2 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 reference_window_stop=0 logging_level=INFO log_to_file=null help=false version=false variant=[(RodBindingCollection [(RodBinding name=variant source=1.gvcf), (RodBinding name=variant2 source=2.gvcf), (RodBinding name=variant3 source=3.gvcf), (RodBinding name=variant4 source=4.gvcf)])] out=Project.test_new_1.vcf includeNonVariantSites=false uniquifySamples=false annotateNDA=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=10 input_prior=[] sample_ploidy=2 annotation=[] group=[Standard] dbsnp=(RodBinding name=dbsnp source=/Reference/dbsnp_137.hg19.vcf) filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
---
##GATKCommandLine.GenotypeGVCFs=<ID=GenotypeGVCFs,Version=3.5-0-g36282e4,Date="Thu Apr 21 11:55:49 BST 2016",Epoch=1461236149722,CommandLineOptions="analysis_type=GenotypeGVCFs input_file=[] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=[/Reference/SeqCap_44Mb_Exome.bed] excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/Reference/ucsc.hg19.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 static_quantized_quals=null round_down_quantized=false disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=true no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=2 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 reference_window_stop=0 logging_level=INFO log_to_file=null help=false version=false variant=[(RodBindingCollection [(RodBinding name=variant source=1.gvcf), (RodBinding name=variant2 source=2.gvcf), (RodBinding name=variant3 source=3.gvcf), (RodBinding name=variant4 source=4.gvcf)])] out=Project.test_new_2.vcf includeNonVariantSites=false uniquifySamples=false annotateNDA=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=10 input_prior=[] sample_ploidy=2 annotation=[] group=[Standard] dbsnp=(RodBinding name=dbsnp source=/Reference/dbsnp_137.hg19.vcf) filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
145,148c145,148
chr1  888639  rs3748596       T       C       5306.35 .       AC=8;AF=1.00;AN=8;DB;DP=150;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=50.20;MQ0=0;QD=29.09;SOR=1.061      GT:AD:DP:GQ:PL  1/1:0,41:41:99:1498,123,0       1/1:0,49:49:99:1714,147,0       1/1:0,26:26:78:930,78,0 1/1:0,34:34:99:1190,102,0
chr1  888659  rs3748597       T       C       4601.35 .       AC=8;AF=1.00;AN=8;DB;DP=132;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=52.14;MQ0=0;QD=32.93;SOR=1.420      GT:AD:DP:GQ:PL  1/1:0,32:32:96:1173,96,0        1/1:0,43:43:99:1522,129,0       1/1:0,28:28:84:896,84,0 1/1:0,28:28:84:1036,84,0
chr1  889158  rs56262069      G       C       4026.35 .       AC=8;AF=1.00;AN=8;DB;DP=89;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=57.77;MQ0=0;QD=30.83;SOR=1.408       GT:AD:DP:GQ:PGT:PID:PL  1/1:0,22:22:69:1|1:889158_G_C:1024,69,0 1/1:0,26:26:78:1|1:889158_G_C:1170,78,0 1/1:0,18:18:57:1|1:889158_G_C:823,57,0  1/1:0,23:23:69:1|1:889158_G_C:1035,69,0
chr1  889159  rs13302945      A       C       4026.35 .       AC=8;AF=1.00;AN=8;DB;DP=91;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=57.77;MQ0=0;QD=29.67;SOR=1.389       GT:AD:DP:GQ:PGT:PID:PL  1/1:0,23:23:69:1|1:889158_G_C:1024,69,0 1/1:0,26:26:78:1|1:889158_G_C:1170,78,0 1/1:0,19:19:57:1|1:889158_G_C:823,57,0  1/1:0,23:23:69:1|1:889158_G_C:1035,69,0
---
chr1  888639  rs3748596       T       C       5306.35 .       AC=8;AF=1.00;AN=8;DB;DP=150;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=50.20;MQ0=0;QD=30.83;SOR=1.061      GT:AD:DP:GQ:PL  1/1:0,41:41:99:1498,123,0       1/1:0,49:49:99:1714,147,0       1/1:0,26:26:78:930,78,0 1/1:0,34:34:99:1190,102,0
chr1  888659  rs3748597       T       C       4601.35 .       AC=8;AF=1.00;AN=8;DB;DP=132;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=52.14;MQ0=0;QD=29.67;SOR=1.420      GT:AD:DP:GQ:PL  1/1:0,32:32:96:1173,96,0        1/1:0,43:43:99:1522,129,0       1/1:0,28:28:84:896,84,0 1/1:0,28:28:84:1036,84,0
chr1  889158  rs56262069      G       C       4026.35 .       AC=8;AF=1.00;AN=8;DB;DP=89;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=57.77;MQ0=0;QD=29.09;SOR=1.408       GT:AD:DP:GQ:PGT:PID:PL  1/1:0,22:22:69:1|1:889158_G_C:1024,69,0 1/1:0,26:26:78:1|1:889158_G_C:1170,78,0 1/1:0,18:18:57:1|1:889158_G_C:823,57,0  1/1:0,23:23:69:1|1:889158_G_C:1035,69,0
chr1  889159  rs13302945      A       C       4026.35 .       AC=8;AF=1.00;AN=8;DB;DP=91;ExcessHet=3.0103;FS=0.000;MLEAC=8;MLEAF=1.00;MQ=57.77;MQ0=0;QD=32.93;SOR=1.389       GT:AD:DP:GQ:PGT:PID:PL  1/1:0,23:23:69:1|1:889158_G_C:1024,69,0 1/1:0,26:26:78:1|1:889158_G_C:1170,78,0 1/1:0,19:19:57:1|1:889158_G_C:823,57,0  1/1:0,23:23:69:1|1:889158_G_C:1035,69,0

Thanks in forward!

Best Answer

Answers

  • JverlouwJverlouw Erasmus MC, RotterdamMember

    Hello Geraldine.

    This does explain it, sort of.

    However, I do notice this effect as well when performing the GenotypesGVCFs on a single thread. This still has to do with the downsampling I assume?

    The difference is not much thruth be told, just around ~12 SNVs that are passing VQSR in just one of the two tests, on a total set of 50k SNVs. Still, such variation in the datasets makes me curious, so that's why I asked.

    Thanks for the explanation!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Yep, it's a fair question. We need to make this into a FAQ article at some point.

Sign In or Register to comment.