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.

Understanding GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0 and GT:DP:GQ:MIN_DP:PL 0/0:1:0:0:0,0,0

mmokrejsmmokrejs Czech RepublicMember

Hi GATK team,
I used the following command to create a per-sample GVCF:

java -jar GenomeAnalysisTK-3.7/GenomeAnalysisTK.jar -T HaplotypeCaller --genotyping_mode DISCOVERY --emitRefConfidence GVCF --pcr_indel_model CONSERVATIVE --sample_name
mysample -R hs38DH.fa -L S04380110__Agilent_SureSelect_Human_All_Exon_V5/hg38/S04380110_Covered.bed --interval_padding 100 -I mysample.bam --dbsnp ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/GATK/00-All.vcf.gz -o mysample.HaplotypeCaller.g.vcf

chr1    980363  .       C       <NON_REF>       .       .       END=980375      GT:DP:GQ:MIN_DP:PL      0/0:4:9:3:0,9,107
chr1    980376  .       C       <NON_REF>       .       .       END=980385      GT:DP:GQ:MIN_DP:PL      0/0:3:6:3:0,6,90
chr1    980386  .       T       <NON_REF>       .       .       END=980394      GT:DP:GQ:MIN_DP:PL      0/0:2:3:2:0,3,45
chr1    980395  .       G       <NON_REF>       .       .       END=980421      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
chr1    980422  .       T       <NON_REF>       .       .       END=980448      GT:DP:GQ:MIN_DP:PL      0/0:1:3:1:0,3,22
chr1    980449  .       C       <NON_REF>       .       .       END=980472      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,47
chr1    980473  .       C       A,<NON_REF>     22.58   .       DP=2;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=7200.00  GT:AD:DP:GQ:PL:SB       1/1:0,2,0:2:6:49,6,0,49,6,49:0,0,1,1
chr1    980474  .       G       <NON_REF>       .       .       END=980485      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,49
chr1    980486  .       G       <NON_REF>       .       .       END=980511      GT:DP:GQ:MIN_DP:PL      0/0:1:3:1:0,3,29
chr1    980512  .       G       <NON_REF>       .       .       END=980531      GT:DP:GQ:MIN_DP:PL      0/0:1:0:0:0,0,0
chr1    980532  .       G       <NON_REF>       .       .       END=980554      GT:DP:GQ:MIN_DP:PL      0/0:1:3:1:0,3,22
chr1    980555  .       A       <NON_REF>       .       .       END=980576      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,39
chr1    980577  .       A       <NON_REF>       .       .       END=980577      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
chr1    980578  .       C       <NON_REF>       .       .       END=980592      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,47
chr1    980593  .       G       <NON_REF>       .       .       END=980612      GT:DP:GQ:MIN_DP:PL      0/0:2:3:1:0,3,29
chr1    980613  .       C       <NON_REF>       .       .       END=980613      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
chr1    980614  .       T       <NON_REF>       .       .       END=980618      GT:DP:GQ:MIN_DP:PL      0/0:1:3:1:0,3,27
chr1    980619  .       C       <NON_REF>       .       .       END=980681      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,47
chr1    980682  .       G       <NON_REF>       .       .       END=980829      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
chr1    981011  .       T       <NON_REF>       .       .       END=981101      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0

I used the BED file with capture tagrgets to obtain only "gene-coding" region in the GVCF, the sequence outside these regions should yield "no call" information, so "./." in the final VCF file. My next step will be GenotypeGVCF, merging 35 such exomes together. I hope that if all my samples some range will be "GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0" that also teh merged VCF file will encode the whole block as "zeroes". But, you advise against using hard DP filters while GVCF is created.

However, how can I get omitted from the per-sample GVCFs lines like:
chr1 980512 . G . . END=980531 GT:DP:GQ:MIN_DP:PL 0/0:1:0:0:0,0,0
chr1 980532 . G . . END=980554 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,22

Do they have any information value?

What "GT:DP:GQ:MIN_DP:PL 0/0:1:0:0:0,0,0" really means?

I attach a figure of the region (sorry, not the BAM with -bamout option used so it may not be the real BAM file content used by HaplotypeCaller) but I do not believe chr1:980523-980532 gained magically a coverage one out of a sudden. Well, HaplotypeCaller actually pokes into soft- or even hard-clipped portions of reads, so it could be teh case ... but no, all four reads shown in the image snapshot had cigar string, 74M.

All in all, I would prefer much more discarding coverage with DP <7 as those are just mis-mapped reads from bwa before getting the per-sample GVCF out. Is that doable? Will "PrintReads -L file.bed" chop the reads and shorten cigar string, or discard reads at the border completely? I could have fed "bwa mem" with the ranges in BED file as well, now it is a bit late.

Finally, would comment what GenotypeGVCFs is likely to do if some sample GVCFs contain same region with roughly same values while say just one out 35 has a bit higher coverage in this particular region? What will happen?

I attached PDF from a VQSR attempt on this dataset to http://gatkforums.broadinstitute.org/gatk/discussion/9394/should-i-provide-the-exome-target-list-l-argu-even-while-calling-gvcf-file-using-haplotypecaller and I bet you will agree that is a bad result. Therefore, I went to re-do the GVCFs and provide "-L range.bed" to restrict the calling region. But as you can see in this current thread, there are still genome-encoding regions with way too low coverage or even no coverage.

Thank you for your comments.

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mmokrejs
    Hi,

    This should not happen in the latest version. If the GQ is 0, the genotype should be ./. to indicate the tool is not confident in any particular genotype. Can you confirm those sites are not set to ./. in the output of GenotypeGVCFs?

    Thanks,
    Sheila

  • mmokrejsmmokrejs Czech RepublicMember
    ##fileformat=VCFv4.2
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##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=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
    ##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.7-0-gcfedb67,Date="Wed Apr 26 11:21:38 CEST 2017",Epoch=1493198498014,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=ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/hs38DH.fa 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 secondsBetweenProgressUpdates=10 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=16 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=./mysample-10_10-PB/realignedBAM/mysample-10_10-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant2 source=./mysample-18_18-PB/realignedBAM/mysample-18_18-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant3 source=./mysample-19_19-PB/realignedBAM/mysample-19_19-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant4 source=./mysample-20_20-PB/realignedBAM/mysample-20_20-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant5 source=./mysample-21_21-PB/realignedBAM/mysample-21_21-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant6 source=./mysample-23_23-PB/realignedBAM/mysample-23_23-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant7 source=./mysample-33_33-PB/realignedBAM/mysample-33_33-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant8 source=./mysample-37_37-PB/realignedBAM/mysample-37_37-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant9 source=./mysample-38_38-PB/realignedBAM/mysample-38_38-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant10 source=./mysample-40_40-PB/realignedBAM/mysample-40_40-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant11 source=./mysample-41_41-PB/realignedBAM/mysample-41_41-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant12 source=./mysample-42_42-PB/realignedBAM/mysample-42_42-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant13 source=./mysample-50_50-PB/realignedBAM/mysample-50_50-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant14 source=./mysample-51_51-PB/realignedBAM/mysample-51_51-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant15 source=./mysample-54_54-PB/realignedBAM/mysample-54_54-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant16 source=./mysample-59_59-PB/realignedBAM/mysample-59_59-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant17 source=./mysample-61_61-PB/realignedBAM/mysample-61_61-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant18 source=./mysample-62_62-PB/realignedBAM/mysample-62_62-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant19 source=./mysample-63_63-PB/realignedBAM/mysample-63_63-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant20 source=./mysample-64_64-PB_excluded/realignedBAM/mysample-64_64-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant21 source=./mysample-65_65-PB/realignedBAM/mysample-65_65-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant22 source=./mysample-66_66-PB/realignedBAM/mysample-66_66-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant23 source=./mysample-68_68-PB_excluded/realignedBAM/mysample-68_68-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant24 source=./mysample-69_69-PB/realignedBAM/mysample-69_69-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant25 source=./mysample-70_70-PB/realignedBAM/mysample-70_70-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant26 source=./mysample-71_71-PB/realignedBAM/mysample-71_71-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant27 source=./mysample-72_72-PB/realignedBAM/mysample-72_72-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant28 source=./mysample-73_73-PB/realignedBAM/mysample-73_73-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant29 source=./mysample-75_75-PB/realignedBAM/mysample-75_75-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant30 source=./mysample-80_80-PB/realignedBAM/mysample-80_80-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant31 source=./mysample-85_85-PB/realignedBAM/mysample-85_85-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant32 source=./mysample-86_86-PB/realignedBAM/mysample-86_86-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant33 source=./mysample-87_87-PB/realignedBAM/mysample-87_87-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant34 source=./mysample-88_88-PB/realignedBAM/mysample-88_88-PB.bwa.gatk.HaplotypeCaller.g.vcf)]), (RodBindingCollection [(RodBinding name=variant35 source=./mysample-89_89-PB/realignedBAM/mysample-89_89-PB.bwa.gatk.HaplotypeCaller.g.vcf)])] out=MGUS_WES_data_n35/mysample.raw.vcf includeNonVariantSites=false uniquifySamples=false annotateNDA=false useNewAFCalculator=true heterozygosity=0.001 indel_heterozygosity=1.25E-4 heterozygosity_stdev=0.01 standard_min_confidence_threshold_for_calling=10.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 max_genotype_count=1024 max_num_PL_values=100 input_prior=[] sample_ploidy=2 annotation=[] group=[StandardAnnotation] dbsnp=(RodBinding name=dbsnp source=ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/GATK/00-All.vcf.gz) filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
    ##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=3.7-0-gcfedb67,Date="Thu Apr 20 19:46:16 CEST 2017",Epoch=1492710376506,CommandLineOptions="analysis_type=HaplotypeCaller input_file=[mysample-18_18-PB.bwa.postalt.removed_duplicates.realignedtogether.BQSR.namesorted.fixmate.calmd.sorted.bam] showFullBamList=false read_buffer_size=null read_filter=[] disable_read_filter=[] intervals=[S04380110__Agilent_SureSelect_Human_All_Exon_V5/hg38/S04380110_Covered.bed] excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=100 reference_sequence=ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/hs38DH.fa 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 secondsBetweenProgressUpdates=10 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 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=dbsnp source=ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/GATK/00-All.vcf.gz) 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 useNewAFCalculator=false heterozygosity=0.001 indel_heterozygosity=1.25E-4 heterozygosity_stdev=0.01 standard_min_confidence_threshold_for_calling=-0.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 max_genotype_count=1024 max_num_PL_values=100 input_prior=[] sample_ploidy=2 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=false 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 maxReadsInMemoryPerSample=30000 maxTotalReadsInMemory=10000000 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">
    
    
    chr1    980565  .       C       A       35.45   .       AC=2;AF=0.053;AN=38;DP=43;ExcessHet=0.0591;FS=0.000;InbreedingCoeff=0.1767;MLEAC=2;MLEAF=0.053;MQ=60.00;QD=17.73;SOR=0.693      GT:AD:DP:GQ:PL  0/0:2,0:2:6:0,6,39      1/1:0,2:2:6:48,6,0      0/0:2,0:2:6:0,6,57      ./.:0,0:0:.:0,0,0       0/0:3,0:3:9:0,9,88
          ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,49      ./.:0,0:0:.:0,0,0       0/0:1,0:1:3:0,3,21      0/0:2,0:2:6:0,6,37      0/0:1,0:1:3:0,3,32      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,61      0/0:4,0:4:9:0,9,135     0/0:4,0:4:12:0,12,99    ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       
    ./.:0,0:0:.:0,0,0       0/0:5,0:5:15:0,15,113   0/0:1,0:1:3:0,3,28      0/0:1,0:1:3:0,3,30      0/0:4,0:4:12:0,12,89    ./.:0,0:0:.:0,0,0       0/0:3,0:3:9:0,9,72      0/0:1,0:1:3:0,3,19      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,53      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       ./.:0,
    0:0:.:0,0,0       0/0:1,0:1:3:0,3,23      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0
    

    In mysamples.recalibrated_variants.vcf I have:

    chr1    980565  .       C       A       35.45   PASS    AC=2;AF=0.053;AN=38;DP=43;ExcessHet=0.0591;FS=0.000;InbreedingCoeff=0.1767;MLEAC=2;MLEAF=0.053;MQ=60.00;QD=17.73;SOR=0.693;VQSLOD=5.74;culprit=SOR      GT:AD:DP:GQ:PL  0/0:2,0:2:6:0,6,39      1/1:0,2:2:6:48,6,0      0/0:2,0:2:6:0,6,57      ./.:0,0:0:.:0,0,0
           0/0:3,0:3:9:0,9,88      ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,49      ./.:0,0:0:.:0,0,0       0/0:1,0:1:3:0,3,21      0/0:2,0:2:6:0,6,37      0/0:1,0:1:3:0,3,32      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,61      0/0:4,0:4:9:0,9,135     0/0:4,0:4:12:0,12,99    ./.:0,0:0:.:0,0,0
           ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:5,0:5:15:0,15,113   0/0:1,0:1:3:0,3,28      0/0:1,0:1:3:0,3,30      0/0:4,0:4:12:0,12,89    ./.:0,0:0:.:0,0,0       0/0:3,0:3:9:0,9,72      0/0:1,0:1:3:0,3,19      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:2,0:2:6:0,6,53      ./.:0,0:0:.:0,0,0
           ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0       0/0:1,0:1:3:0,3,23      ./.:0,0:0:.:0,0,0       ./.:0,0:0:.:0,0,0
    

    The BED file contains at about this locus:

    chr1    978880  979933  ref|PERM1,ref|NM_001291367,ref|NM_001291366,ens|ENST00000341290,ens|ENST00000433179,mRNA|AK123855,mRNA|KF150175
    chr1    980108  980358  ref|PERM1,ref|NM_001291367,ref|NM_001291366,ens|ENST00000341290,ens|ENST00000433179,mRNA|AK123855,mRNA|KF150175
    chr1    980362  980729  ref|PERM1,ref|NM_001291367,ref|NM_001291366,ens|ENST00000341290,ens|ENST00000433179,mRNA|AK123855,mRNA|KF150175
    chr1    981110  981283  ref|PERM1,ref|NM_001291367,ref|NM_001291366,ens|ENST00000341290,ens|ENST00000433179,mRNA|AK123855,mRNA|KF150175
    

    So, this is really low-covered region either due to sequence capture or due to just low representation in the sequencing library. However, does this impose any problem to the VQSR approach? Note, the PDF files attached to http://gatk.vanillaforums.com/discussion/comment/38257 showing the filtration went badly. Initially I had much more such lines in the VCF file, now thanks to -L file.bed I have much less of these but the VQSR results seem still just bad.

    Let's concern here what the GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0 and GT:DP:GQ:MIN_DP:PL 0/0:1:0:0:0,0,0 say and whether it would be more meaningful to omit such rows from the VCF output. What is wrong with my VQSR attempt I would better discuss at http://gatk.vanillaforums.com/discussion/comment/38257 . Thank you.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mmokrejs
    Hi,

    The low-covered regions should not pose an problem to VQSR. Those are usually at the edge of the capture regions. VariantRecalibrator does not take into account the genotype-level information. It only uses the site-level information. Most site-level annotations are calculated based on the samples that actually have a variant at the site, so the no-calls and hom-ref genotypes will not have any impact.

    As for the 0/0 genotypes with 0 GQ, I don't see any of those in the last post. The GQ 0 sites are indeed set to no-calls.

    What do you mean you had many lines in the VCF but now have much less? What is the BED file you are using?

    Thanks,
    Sheila

  • mmokrejsmmokrejs Czech RepublicMember

    @Sheila

    I see, I showed the merged multi-sample VCF files. Here are the per-sample GVCFs:

    $ find . -name \*-PB.bwa.gatk.HaplotypeCaller.g.vcf | while read f; do grep -H "^chr1" $f | awk '{if ($2==980395) print}'; done
    ./mysample-54_54-PB/realignedBAM/mysample-54_54-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980494      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
    ./mysample-18_18-PB/realignedBAM/mysample-18_18-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980423      GT:DP:GQ:MIN_DP:PL      0/0:2:6:2:0,6,47
    ./mysample-10_10-PB/realignedBAM/mysample-10_10-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980421      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
    ./mysample-71_71-PB/realignedBAM/mysample-71_71-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980396      GT:DP:GQ:MIN_DP:PL      0/0:2:3:2:0,3,45
    ./mysample-61_61-PB/realignedBAM/mysample-61_61-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980565      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
    ./mysample-87_87-PB/realignedBAM/mysample-87_87-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980397      GT:DP:GQ:MIN_DP:PL      0/0:6:6:5:0,6,90
    ./mysample-73_73-PB/realignedBAM/mysample-73_73-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980398      GT:DP:GQ:MIN_DP:PL      0/0:4:6:4:0,6,90
    ./mysample-89_89-PB/realignedBAM/mysample-89_89-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980397      GT:DP:GQ:MIN_DP:PL      0/0:6:15:6:0,15,225
    ./mysample-51_51-PB/realignedBAM/mysample-51_51-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980399      GT:DP:GQ:MIN_DP:PL      0/0:3:6:3:0,6,90
    ./mysample-86_86-PB/realignedBAM/mysample-86_86-PB.bwa.gatk.HaplotypeCaller.g.vcf:chr1    980395  .       G       <NON_REF>       .       .       END=980829      GT:DP:GQ:MIN_DP:PL      0/0:0:0:0:0,0,0
    

    In other samples the chr1:980395 is buried the GVCF block so I cannot easily look for this position.

    Initially I created g.vcf files from files without using -L option to restrict regions to Agilent Human Exome v5 ranges. That is why I had so many 0/0 genotypes. After redoing the HaplotypeCaller step with a BED file the number of lines with practically no coverage shrunk down. Those remaining are occasionally inside of captured regions, something I have to live with. So those 0/0:0:0:0:0,0,0 were coming from the non-targeted regions. However, some are still left and as you say there should be ./. instead.

  • mmokrejsmmokrejs Czech RepublicMember

    Fine, so how about my original question: However, how can I get omitted from the per-sample GVCFs lines like: ...? Can I get such lines not even included in the per-sample GVCF? For size reasons at least. The other question whether such GQ0 hom-ref fool downstream analyses or not seems best answred with "Should not provided such samples get represented in the final VCF with no calls: ./.".

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    The only way to avoid having these lines in the GVCF is to exclude those regions from analysis in the first place. Eg if you do a coverage analysis first and find which regions aren't callable anyway across all or enough of your samples to warrant exclusion, you can specify those regions to be excluded.
  • mmokrejsmmokrejs Czech RepublicMember
Sign In or Register to comment.