Bug Bulletin: The GenomeLocPArser error in SplitNCigarReads has been fixed; if you encounter it, use the latest nightly build.

Unified Genotyper Allele Depth (AD) filed for samples

hanshans Posts: 3Member

Hello,

I found some strange entries for indels in my VCF file created by the Unified Genotyper. For example:

4 184513470 . TC T 4009 PASS AC=4;AF=0.250;AN=16;BaseQRankSum=1.972;DP=1315;DS;FS=3.466;HaplotypeScore=537.6937;MLEAC=4;MLEAF=0.250;MQ=52.55;MQ0=0;MQRankSum=-10.581;QD=4.55;ReadPosRankSum=-10.128;SB=-3.500e+01;set=variant2 GT:AD:DP:GQ:PL 0/1:230,0:239:99:282,0,5011 0/0:92,0:95:99:0,133,2435

The first sample has genotype 0/1 with a good GQ value. However, according the allele depth field, there is no read supporting the deletion. When I look at the reads using the IGV, I find some reads supporting the deletion for the first sample (and even some for the second one).

Moreover, when I looked at the AD values for SNPs, I noticed the the sum of all AD values is much less than the coverage shown in the IGV. I filtered duplicated reads in the IGV.

Can someone please give an explanation? This link http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_annotator_DepthPerAlleleBySample.html explains the difference between AD and DP, but does not help in my case.

Best greetings, Hans-Ulrich

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,176Administrator, GATK Developer admin

    Hi Hans, what version of the GATK are you using?

    Geraldine Van der Auwera, PhD

  • hanshans Posts: 3Member

    Sorry, I forgot to give the version number. It is "The Genome Analysis Toolkit (GATK) v2.1-8-g5efb575, Compiled 2012/08/30 14:22:17"

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,176Administrator, GATK Developer admin

    Ah, you'll need to upgrade to the latest version. The issue you're having has been fixed.

    Geraldine Van der Auwera, PhD

  • hanshans Posts: 3Member

    OK, I will try the new version. Thanks!

  • dolsemtldolsemtl Posts: 1

    Hi, I got the same entries with GATK-lite. Did you solve the problem with the latest version. If so, which version you did. Thanks.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,176Administrator, GATK Developer admin

    Hi @dolsemtl, I recommend you use the very latest version (currently 2.3-9).

    Geraldine Van der Auwera, PhD

  • alimanaliman Oxford, UKPosts: 4Member

    Hi, using GATK 2.3-9 I'm finding that the allele depths (AD field in calldata) for deletions called with the UnifiedGenotyper (ploidy=1) are zero or very low for both ref and alt where there are clearly many reads supporting one or the other allele. Here's an example:

    Pf3D7_04_v3 99839 . CTATATATATATATATA C 55455.68 PASS AC=13;AF=0.591;AN=22;BaseCounts=75,1255,0,1;BaseQRankSum=-3.665;DP=1343;FS=0.0;HRun=0;HaplotypeScore=65.1358;IndelType=DEL.NOVEL_10orMore.;LowMQ=0.0067,0.0067,1343.0;MLEAC=13;MLEAF=0.591;MQ=53.45;MQ0=0;MQ0Fraction=0.0067;MQRankSum=-2.757;QD=61.89;RPA=19,11;RU=TA;ReadPosRankSum=1.839;SB=-23010.0;STR;RegionType=Conserved;DPSoftClipped=130;DPMateSameStrand=0;DPFaceAway=0;DPAll=1158;DPProperPair=1153;DPMateOtherChrom=2;DPDuplicate=0;DPMateUnmapped=3;DPNormedByGC=0.871987951807;DPNormed=0.535120147874;DPPercentile=15;DPPercentileByGC=38;StdTlen=299;MeanTlen=26;RMSTlen=300;SoftClippingBias=7;UQ=45;GC100=7.0;GC300=11.0;GC500=14.0;HProx5=24;HProx10=275;HProx15=485;DPProperPairFraction=0.9957;DPMateUnmappedFraction=0.0026;DPMateOtherChromFraction=0.0017;DPMateSameStrandFraction=0.0;DPFaceAwayFraction=0.0;DPSoftClippedFraction=0.1123;DPDuplicateFraction=0.0;DPBiasStudent=9;DPBiasWelch=9;DPBiasRankSum=8 GT:AD:DP:GQ:MLPSAC:MLPSAF:MQ0:PL 1:0,1:62:99:1:1.0:0:3233,0 0:0,0:53:99:0:0.0:0:0,4619 0:0,0:75:99:0:0.0:0:0,6389 0:0,0:51:99:0:0.0:0:0,3802 1:0,0:85:99:1:1.0:0:6233,0 0:0,0:82:99:0:0.0:0:0,8012 1:0,2:117:99:1:1.0:1:6877,0 0:0,0:12:99:0:0.0:0:0,586 1:0,2:67:99:1:1.0:0:4769,0 1:0,1:34:99:1:1.0:0:2814,0 0:0,0:56:99:0:0.0:1:0,5989 0:0,0:15:99:0:0.0:0:0,1766 1:0,1:62:99:1:1.0:1:3643,0 1:0,2:50:99:1:1.0:2:2681,0 1:0,4:79:99:1:1.0:1:5223,0 0:0,0:33:99:0:0.0:0:0,2154 0:0,0:65:99:0:0.0:1:0,6721 1:0,0:25:99:1:1.0:0:2167,0 1:0,1:57:99:1:1.0:0:4339,0 1:0,1:116:99:1:1.0:2:7482,0 1:0,4:77:99:1:1.0:0:4884,0 1:0,1:49:99:1:1.0:0:1201,0

    Here's the VCF header:

    UnifiedGenotyper="analysis_type=UnifiedGenotyper input_file=[./data/3d7_v3/bwa_n0.01_k4_l32/alignment_realigned_merged/3d7_hb3.bam] read_buffer_size=null phone_home=STANDARD gatk_key=null tag=NA read_filter=[] intervals=[Pf3D7_04_v3] excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=./data/genome/sanger/version3/September_2012/Pf3D7_v3.fa nonDeterministicRandomSeed=false disableRandomization=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=NONE downsample_to_fraction=null downsample_to_coverage=null use_legacy_downsampler=false baq=OFF baqGapOpenPenalty=40.0 fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false performanceLog=null useOriginalQualities=false BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 defaultBaseQualities=-1 validation_strictness=SILENT remove_program_records=false keep_program_records=false unsafe=null 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 logging_level=INFO log_to_file=null help=false genotype_likelihoods_model=INDEL pcr_error_rate=1.0E-4 computeSLOD=true annotateNDA=false pair_hmm_implementation=ORIGINAL min_base_quality_score=17 max_deletion_fraction=0.05 min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indel_heterozygosity=1.25E-4 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null sample_ploidy=1 min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 min_reference_depth=100 exclude_filtered_reference_sites=false heterozygosity=0.001 genotyping_mode=DISCOVERY output_mode=EMIT_VARIANTS_ONLY standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 alleles=(RodBinding name= source=UNBOUND) max_alternate_alleles=6 p_nonref_model=EXACT_GENERAL_PLOIDY contamination_fraction_to_filter=0.0 logRemovedReadsFromContaminationFiltering=null exactcallslog=null dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub no_cmdline_in_header=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub bcf=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub debug_file=null metrics_file=null annotation=[IndelType, BaseCounts, BaseQualityRankSumTest, ChromosomeCounts, ClippingRankSumTest, DepthOfCoverage, DepthPerAlleleBySample, FisherStrand, HaplotypeScore, HomopolymerRun, LowMQ, MappingQualityRankSumTest, MappingQualityZero, MappingQualityZeroBySample, MappingQualityZeroFraction, QualByDepth, RMSMappingQuality, ReadPosRankSumTest, SampleList, SpanningDeletions] excludeAnnotation=[] filter_mismatching_base_and_quals=false"

    Cheers,

    Alistair

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,176Administrator, GATK Developer admin

    Hi Alistair,

    There's not much we can say about this without seeing some pileups or other data. There's a lot of things that could explain why the supporting reads you see aren't being counted. I recommend you check out this article:

    http://www.broadinstitute.org/gatk/guide/article?id=1235

    If you still think the AD numbers are incorrect, take an IGV screenshot of your example region and show us which reads you think aren't being counted.

    Geraldine Van der Auwera, PhD

  • celzingacelzinga Posts: 1Member
    edited March 2013

    Hello,

    I'm seeing a similar problem (using 2.3.9), and I'm noticing a pattern: it seems only to occur when there is a SNP at the same location for one of the samples in the UG run. Oddly enough, it happens regardless of whether I'm using -glm BOTH or -glm INDEL. If I omit the sample(s) that contain the SNP, then the insertion/deletion reads do get counted in AD.

    Here's a part of the VCF from a UG run that included a sample with a SNP at the location (the sample with the SNP is not shown in the snippet below): 6 110053824 . G T 2757.13 . AC=1;AF=0.031;AN=32;BaseQRankSum=-14.015;DP=3471;Dels=0.00;FS=48.970;HaplotypeScore=191.7272;InbreedingCoeff=-0.0360;MLEAC=1;MLEAF=0.031;MQ=56.65;MQ0=2;MQRankSum=3.472;QD=11.03;ReadPosRankSum=-6.390 GT:AD:DP:GQ:PL 0/0:236,14:250:99:0,526,6571 6 110053824 rs150311668 G GT 371313.22 . AC=29;AF=0.906;AN=32;DB;DP=3471;FS=2.505;HaplotypeScore=281.0531;InbreedingCoeff=-0.1056;MLEAC=29;MLEAF=0.906;MQ=56.65;MQ0=0;QD=106.98;RPA=13,14;RU=T;STR GT:AD:DP:GQ:PL 1/1:0,0:250:99:7049,428,0

    Then, when I omitted the sample with a SNP at that location for the UG run, I get the insertion count in AD (using the same UG settings): 6 110053824 rs150311668 G GT 32729.73 . AC=2;AF=1.00;AN=2;BaseQRankSum=-0.234;DB;DP=250;FS=4.614;HaplotypeScore=325.8291;MLEAC=2;MLEAF=1.00;MQ=57.69;MQ0=0;MQRankSum=1.010;QD=130.92;RPA=13,14;RU=T;ReadPosRankSum=-0.688;STR GT:AD:DP:GQ:PL 1/1:10,212:250:99:6965,379,0

    I had a sample with both the SNP and the insertion at that location, and that one (surprisingly) did count the insertions in AD. Probably a different issue (due in part, I'm guessing, to the SNP and indel callers being different processes), but the 0/1 call for the SNP and 1/1 call for the insertion appear to be in conflict given it's a single position: 6 110053824 . G T 2995.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-2.852;DP=250;Dels=0.02;FS=3.677;HaplotypeScore=221.5300;MLEAC=1;MLEAF=0.500;MQ=55.11;MQ0=1;MQRankSum=-1.532;QD=11.98;ReadPosRankSum=-2.522 GT:AD:DP:GQ:PL 0/1:125,121:246:99:3024,0,3132 6 110053824 rs150311668 G GT 5077.73 . AC=2;AF=1.00;AN=2;BaseQRankSum=-3.673;DB;DP=250;FS=6.039;HaplotypeScore=789.8322;MLEAC=2;MLEAF=1.00;MQ=55.11;MQ0=0;MQRankSum=0.347;QD=20.31;RPA=13,14;RU=T;ReadPosRankSum=-0.080;STR GT:AD:DP:GQ:PL 1/1:35,187:249:92:5115,92,0

    all the best, Chris

    Post edited by celzinga on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,176Administrator, GATK Developer admin

    Can you try with 2.4 and tell us if the problem persists?

    Geraldine Van der Auwera, PhD

Sign In or Register to comment.