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.

Bug in output VariantFiltration with --genotypeFilterExpression

Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭
in Ask the GATK team

When using VariantFilteration with genotype filters, the FORMAT line for the genotype filter is wrong:
##FORMAT=<ID=FT,Number=1,Type=String,Description="Genotype-level filter">

It should have "Number=." like so:
##FORMAT=<ID=FT,Number=.,Type=String,Description="Genotype-level filter">

Since any genotype can have 0 or more filters associated with it. From the VCFv4.2 doc: "If the number of possible values varies, is unknown, or is unbounded, then this value should be ‘.’"

Its a small issue, I noticed it when pyvcf didn't properly parse the genotype filters to a list. On the plus side, it should be a one character bugfix, which is always nice.

To replicate, I used this command to trigger two filters on one genotype. This occurs both in gatk3.6 and gatk3.5.

gatk -T VariantFiltration \
-R $ref \
--variant $vcf_file \
--filterExpression "QUAL<1000.0" \
--filterName "Q1000" \
--genotypeFilterExpression "DP<20" \
--genotypeFilterName "DP20" \
--genotypeFilterExpression "DP<20" \
--genotypeFilterName "DP20_v2" \
-o $output \
--setFilteredGtToNocall

Tagged:
· Share on Twitter Share on Google+

Best Answer

Answers

  • Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭

    Awesome, thanks a lot!

    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @Redmar_van_den_Berg, this has been fixed in htsjdk and the version of htsjdk was revved in GATK, so the latest nightlies have the fix (since Sept 20 actually).

    · Share on Twitter Share on Google+
  • Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭

    Hi @Geraldine_VdAuwera, thanks for getting back to me.

    I've tested the latest nightly build but I think something went wrong, the FT annotation still shows "Number=1".

    Nightly build:

    $ java -jar ~/Downloads/GenomeAnalysisTK.jar --version
nightly-2016-10-28-gaca5d7b

    Snipped from filtered VCF file, ##GATKCommandLine.VariantFiltration shows the nightly build (I used a vcf file generated with stable gatk).

    ##FORMAT=<ID=FT,Number=1,Type=String,Description="Genotype-level filter">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=RGQ,Number=1,Type=Integer,Description="Unconditional reference genotype confidence, encoded as a phred quality -10*log10 p(genotype call is wrong)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##GATKCommandLine.GenotypeGVCFs=<ID=GenotypeGVCFs,Version=3.6-0-g89b7209,Date="Wed Oct 26 23:17:42 CEST 2016",Epoch=1477516662230,CommandLineOptions="analysis_type=GenotypeGVCFs input_file=[] showFullBamList=false read_buffer_size=null read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=cj_reference.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=false 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=3 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 phone_home= gatk_key=null tag=NA logging_level=INFO log_to_file=null help=false version=false variant=[(RodBindingCollection [(RodBinding name=variant source=alignments/cj10_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant2 source=alignments/cj11_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant3 source=alignments/cj12_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant4 source=alignments/cj14_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant5 source=alignments/cj15_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant6 source=alignments/cj16_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant7 source=alignments/cj1_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant8 source=alignments/cj2_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant9 source=alignments/cj3_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant10 source=alignments/cj4_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant11 source=alignments/cj5_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant12 source=alignments/cj6_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant13 source=alignments/cj7_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant14 source=alignments/cj8_sorted_dedup_realn.g.vcf)]), (RodBindingCollection [(RodBinding name=variant15 source=alignments/cj9_sorted_dedup_realn.g.vcf)])] out=/data/GMI 2016/cj/snps/joint_snips.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=6 max_num_PL_values=100 input_prior=[] sample_ploidy=2 annotation=[] group=[StandardAnnotation] dbsnp=(RodBinding name= source=UNBOUND) filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=3.6-0-g89b7209,Date="Wed Oct 26 22:49:09 CEST 2016",Epoch=1477514949365,CommandLineOptions="analysis_type=HaplotypeCaller input_file=[alignments/cj11_sorted_dedup_realn.bam] showFullBamList=false read_buffer_size=null read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=cj_reference.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=500 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=false 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=1 num_cpu_threads_per_data_thread=3 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 phone_home= gatk_key=null tag=NA logging_level=INFO log_to_file=null help=false version=false likelihoodCalculationEngine=PairHMM heterogeneousKmerSizeResolution=COMBO_MIN dbsnp=(RodBinding name= source=UNBOUND) dontTrimActiveRegions=false maxDiscARExtension=25 maxGGAARExtension=300 paddingAroundIndels=150 paddingAroundSNPs=20 comp=[] annotation=[StrandBiasBySample] excludeAnnotation=[ChromosomeCounts, FisherStrand, StrandOddsRatio, QualByDepth] group=[StandardAnnotation, StandardHCAnnotation] debug=false useFilteredReadsForAnnotations=false emitRefConfidence=GVCF bamOutput=null bamWriterType=CALLED_HAPLOTYPES emitDroppedReads=false disableOptimizations=false annotateNDA=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=-0.0 standard_min_confidence_threshold_for_emitting=-0.0 max_alternate_alleles=6 max_num_PL_values=100 input_prior=[] sample_ploidy=1 genotyping_mode=DISCOVERY alleles=(RodBinding name= source=UNBOUND) contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=null exactcallslog=null output_mode=EMIT_VARIANTS_ONLY allSitePLs=true gcpHMM=10 pair_hmm_implementation=VECTOR_LOGLESS_CACHING pair_hmm_sub_implementation=ENABLE_ALL always_load_vector_logless_PairHMM_lib=false phredScaledGlobalReadMismappingRate=45 noFpga=false sample_name=null kmerSize=[10, 25] dontIncreaseKmerSizesForCycles=false allowNonUniqueKmersInRef=false numPruningSamples=1 recoverDanglingHeads=false doNotRecoverDanglingBranches=false minDanglingBranchLength=4 consensus=false maxNumHaplotypesInPopulation=128 errorCorrectKmers=false minPruning=2 debugGraphTransformations=false allowCyclesInKmerGraphToGeneratePaths=false graphOutput=null kmerLengthForReadErrorCorrection=25 minObservationsForKmerToBeSolid=20 GVCFGQBands=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 70, 80, 90, 99] indelSizeToEliminateInRefModel=10 min_base_quality_score=10 includeUmappedReads=false useAllelesTrigger=false doNotRunPhysicalPhasing=true keepRG=null justDetermineActiveRegions=false dontGenotype=false dontUseSoftClippedBases=false captureAssemblyFailureBAM=false errorCorrectReads=false pcr_indel_model=CONSERVATIVE maxReadsInRegionPerSample=10000 minReadsPerAlignmentStart=10 mergeVariantsViaLD=false activityProfileOut=null activeRegionOut=null activeRegionIn=null activeRegionExtension=null forceActive=false activeRegionMaxSize=null bandPassSigma=null maxProbPropagationDistance=50 activeProbabilityThreshold=0.002 min_mapping_quality_score=20 filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
##GATKCommandLine.VariantFiltration=<ID=VariantFiltration,Version=nightly-2016-10-28-gaca5d7b,Date="Fri Oct 28 11:47:03 CEST 2016",Epoch=1477648023585,CommandLineOptions="analysis_type=VariantFiltration input_file=[] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=../cj_reference.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 static_quantized_quals=null round_down_quantized=false quantize_quals=0 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=false 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=1 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=(RodBinding name=variant source=joint_snips.vcf) mask=(RodBinding name= source=UNBOUND) out=/data/GMI 2016/cj/snps/nightly_out filterExpression=[QUAL<1000.0] filterName=[Q1000] genotypeFilterExpression=[DP<20, DP<20] genotypeFilterName=[DP20, DP20_v2] clusterSize=3 clusterWindowSize=0 maskExtension=0 maskName=Mask filterNotInMask=false missingValuesInExpressionsShouldEvaluateAsFailing=false invalidatePreviousFilters=false invertFilterExpression=false invertGenotypeFilterExpression=false setFilteredGtToNocall=true filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    Were you running on a vcf that already had the annotation? GATK won't replace/update a definition line that's already there unless you request it (I believe we added some logic to that effect recently).
    · Share on Twitter Share on Google+
  • Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭

    @Geraldine_VdAuwera

    The unfiltered vcf file did not have an FT annotation. I can see that the nightly build leaves existing modified annotations alone though, while gatk3.6 overwrites them, so that change has landed in the nightly build.

    · Share on Twitter Share on Google+
  • SheilaSheila Broad InstituteMember, Broadie admin

    @Redmar_van_den_Berg
    Hi,

    Hmm. I just tested the latest nightly build, and I get the correct output. However, when I test the latest stable version (3.6), I get the "buggy" output. Can you post the exact command you ran? Also, can you please try the latest nightly build from last night? That is the one I am using and I can see a difference in the outputs.

    Thanks,
    Sheila

    · Share on Twitter Share on Google+
  • Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭

    @Sheila I don't know whats going on...

    Excerpt from .bash_history that shows the exact commands I used

     2007  export ref=../cj_reference.fasta
 2008  export vcf_file=joint_snips.vcf
 2009  export output=nightly.vcf
 2010  vim joint_snips.vcf
 2011  java -jar ~/Downloads/GenomeAnalysisTK.jar --version
 2012  java -jar ~/Downloads/GenomeAnalysisTK.jar -T VariantFiltration  -R $ref  --variant $vcf_file  --filterExpression "QUAL<1000.0"  --filterName "Q1000"  --genotypeFilterExpression "DP<20"  --genotypeFilterName "DP20"  --genotypeFilterExpression "DP<20"  --genotypeFilterName "DP20_v2"  -o $output  --setFilteredGtToNocall
 2013  vim $output

    First, I checked that the $vcf_file did not contain a FT annotation

    Then I show the nightly build:

    $ java -jar ~/Downloads/GenomeAnalysisTK.jar --version
nightly-2016-11-02-gaca5d7b

    Then I run the filter command, and check the $output vcf file using vim:

    ##FORMAT=<ID=FT,Number=1,Type=String,Description="Genotype-level filter">
    · Share on Twitter Share on Google+
  • SheilaSheila Broad InstituteMember, Broadie admin

    @Redmar_van_den_Berg
    Hi,

    Can you please post the before and after VCF header line?

    Thanks
    Sheila

    · Share on Twitter Share on Google+
  • Redmar_van_den_BergRedmar_van_den_Berg Member ✭✭

    @Sheila
    Hi,

    What do you mean with the before and after VCF header line? The FT header gets added after the variant filtration step, so it is not present in the unfiltered VCF file.
    However, when it gets added it shows

    ##FORMAT=<ID=FT,Number=1,Type=String,Description="Genotype-level filter">

    Which should be

    ##FORMAT=<ID=FT,Number=.,Type=String,Description="Genotype-level filter">
    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    We just wanted to be extra sure that the FT line wasn't already in there from a previous operation. Meanwhile I tried this myself on an unfiltered vcf and get the same result you do. There's clearly something unexpected going on so we'll check everything and get back to you.

    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    OK, we've determined that there is a problem with our build system that is causing a stale version to be packaged as a nightly build. We didn't catch this because we do our fix-testing against a post-merge build, which normally should be equivalent to using the latest nightly -- except right now it's not. Engineers are looking into the problem now; in the meantime if you need the fix you can compile from source to get the correct latest build.

    · Share on Twitter Share on Google+
Sign In or Register to comment.