From a783f19ab12060084c9811902365d7629b1631ca Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Wed, 6 Mar 2013 13:45:53 -0500 Subject: [PATCH] Fix for potential HaplotypeCaller bug in annotation ordering -- Annotations were being called on VariantContext that might needed to be trimmed. Simply inverted the order of operations so trimming occurs before the annotations are added. -- Minor cleanup of call to PairHMM in LikelihoodCalculationEngine --- .../walkers/haplotypecaller/GenotypingEngine.java | 13 ++++++++----- .../LikelihoodCalculationEngine.java | 9 ++++++--- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java index 1cfc65581..400de6485 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java @@ -273,16 +273,19 @@ public class GenotypingEngine { final Map alleleReadMap_annotations = ( USE_FILTERED_READ_MAP_FOR_ANNOTATIONS ? alleleReadMap : convertHaplotypeReadMapToAlleleReadMap( haplotypeReadMap, alleleMapper, 0.0, UG_engine.getUAC().contaminationLog ) ); final Map stratifiedReadMap = filterToOnlyOverlappingReads( genomeLocParser, alleleReadMap_annotations, perSampleFilteredReadList, call ); - VariantContext annotatedCall = annotationEngine.annotateContext(stratifiedReadMap, call); + + VariantContext annotatedCall = call; + // TODO -- should be before annotated call, so that QDL works correctly + if( annotatedCall.getAlleles().size() != mergedVC.getAlleles().size() ) { // some alleles were removed so reverseTrimming might be necessary! + annotatedCall = GATKVariantContextUtils.reverseTrimAlleles(annotatedCall); + } + + annotatedCall = annotationEngine.annotateContext(stratifiedReadMap, annotatedCall); // maintain the set of all called haplotypes for ( final Allele calledAllele : call.getAlleles() ) calledHaplotypes.addAll(alleleMapper.get(calledAllele)); - if( annotatedCall.getAlleles().size() != mergedVC.getAlleles().size() ) { // some alleles were removed so reverseTrimming might be necessary! - annotatedCall = GATKVariantContextUtils.reverseTrimAlleles(annotatedCall); - } - returnCalls.add( annotatedCall ); } } diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java index a7d85b969..87b488b3e 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java @@ -151,9 +151,12 @@ public class LikelihoodCalculationEngine { final int haplotypeStart = ( previousHaplotypeSeen == null ? 0 : PairHMM.findFirstPositionWhereHaplotypesDiffer(haplotype.getBases(), previousHaplotypeSeen.getBases()) ); previousHaplotypeSeen = haplotype; - perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype), - pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), read.getReadBases(), - readQuals, readInsQuals, readDelQuals, overallGCP, haplotypeStart, jjj == 0)); + final boolean isFirstHaplotype = jjj == 0; + final double log10l = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), + read.getReadBases(), readQuals, readInsQuals, readDelQuals, + overallGCP, haplotypeStart, isFirstHaplotype); + + perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype), log10l); } } return perReadAlleleLikelihoodMap;