From f78ff08e2b732eddecf3e07214bc42b963856455 Mon Sep 17 00:00:00 2001 From: ebanks Date: Fri, 15 Oct 2010 18:56:45 +0000 Subject: [PATCH] This is less correct than my previous change but it's what UGv1 does and now is not the right time to start mucking with things. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4506 348d0f76-0448-11de-a6fe-93d51630548a --- .../genotyper/AlleleFrequencyCalculationModel.java | 4 +++- .../walkers/genotyper/ExactAFCalculationModel.java | 6 +++--- .../gatk/walkers/genotyper/GridSearchAFEstimation.java | 10 ++++++++-- .../gatk/walkers/genotyper/UnifiedGenotyperEngine.java | 8 +++----- 4 files changed, 17 insertions(+), 11 deletions(-) diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java index 730b9de93..4ce32d971 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java @@ -81,12 +81,14 @@ public abstract class AlleleFrequencyCalculationModel implements Cloneable { * @param GLs genotype likelihoods * @param log10AlleleFrequencyPriors priors * @param log10AlleleFrequencyPosteriors array (pre-allocated) to store results + * @param minFrequencyToCalculate the minimum frequency which needs to be calculated */ public abstract void getLog10PNonRef(RefMetaDataTracker tracker, ReferenceContext ref, Map GLs, double[] log10AlleleFrequencyPriors, - double[] log10AlleleFrequencyPosteriors); + double[] log10AlleleFrequencyPosteriors, + int minFrequencyToCalculate); /** * Can be overridden by concrete subclasses diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/ExactAFCalculationModel.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/ExactAFCalculationModel.java index dfa455a69..c6c8938ee 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/ExactAFCalculationModel.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/ExactAFCalculationModel.java @@ -29,7 +29,6 @@ import org.apache.log4j.Logger; import org.broad.tribble.util.variantcontext.Allele; import org.broad.tribble.util.variantcontext.Genotype; import org.broad.tribble.util.variantcontext.GenotypeLikelihoods; -import org.broad.tribble.util.variantcontext.VariantContext; import org.broad.tribble.vcf.VCFConstants; import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext; import org.broadinstitute.sting.utils.*; @@ -54,7 +53,8 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { ReferenceContext ref, Map GLs, double[] log10AlleleFrequencyPriors, - double[] log10AlleleFrequencyPosteriors) { + double[] log10AlleleFrequencyPosteriors, + int minFrequencyToCalculate) { // Math requires linear math to make efficient updates. double[] alleleFrequencyPriors = MathUtils.normalizeFromLog10(log10AlleleFrequencyPriors); @@ -233,7 +233,7 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { HashMap attributes = new HashMap(); ArrayList myAlleles = new ArrayList(); attributes.put(VCFConstants.DEPTH_KEY, contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size()); - double qual = 0.0; + double qual; double[] posteriors = GLs.get(sample).getPosteriors(); if (bestGTguess == 0) { diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/GridSearchAFEstimation.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/GridSearchAFEstimation.java index edef58d13..00adb3757 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/GridSearchAFEstimation.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/GridSearchAFEstimation.java @@ -48,7 +48,8 @@ public class GridSearchAFEstimation extends AlleleFrequencyCalculationModel { ReferenceContext ref, Map GLs, double[] log10AlleleFrequencyPriors, - double[] log10AlleleFrequencyPosteriors) { + double[] log10AlleleFrequencyPosteriors, + int minFrequencyToCalculate) { initializeAFMatrix(GLs); @@ -68,12 +69,17 @@ public class GridSearchAFEstimation extends AlleleFrequencyCalculationModel { // an optimization to speed up the calculation: if we are beyond the local maximum such // that subsequent likelihoods won't factor into the confidence score, just quit - if ( maxLikelihoodSeen - log10AlleleFrequencyPosteriors[i] > LOG10_OPTIMIZATION_EPSILON ) + if ( i >= minFrequencyToCalculate && maxLikelihoodSeen - log10AlleleFrequencyPosteriors[i] > LOG10_OPTIMIZATION_EPSILON ) return; if ( log10AlleleFrequencyPosteriors[i] > maxLikelihoodSeen ) maxLikelihoodSeen = log10AlleleFrequencyPosteriors[i]; } + + // TODO: we really need to get rid of this and the minFrequencyToCalculate argument + // it's possible that we need to calculate higher frequencies + for (int i = maxAlleleFrequencyToTest; i <= minFrequencyToCalculate; i++) + log10AlleleFrequencyPosteriors[i] = log10AlleleFrequencyPosteriors[maxAlleleFrequencyToTest]; } /** diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java index 57e67764c..55959f7c7 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java @@ -157,7 +157,7 @@ public class UnifiedGenotyperEngine { // 'zero' out the AFs (so that we don't have to worry if not all samples have reads at this position) clearAFarray(log10AlleleFrequencyPosteriors.get()); - afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get()); + afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), 0); // find the most likely frequency int bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get()); @@ -234,9 +234,8 @@ public class UnifiedGenotyperEngine { GLs.clear(); glcm.get().getLikelihoods(tracker, refContext, stratifiedContexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD, genotypePriors, GLs); clearAFarray(log10AlleleFrequencyPosteriors.get()); - afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get()); + afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), bestAFguess); double forwardLog10PofNull = log10AlleleFrequencyPosteriors.get()[0]; - bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get()); double forwardLog10PofF = log10AlleleFrequencyPosteriors.get()[bestAFguess]; //System.out.println("forwardLog10PofNull=" + forwardLog10PofNull + ", forwardLog10PofF=" + forwardLog10PofF); @@ -244,9 +243,8 @@ public class UnifiedGenotyperEngine { GLs.clear(); glcm.get().getLikelihoods(tracker, refContext, stratifiedContexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE, genotypePriors, GLs); clearAFarray(log10AlleleFrequencyPosteriors.get()); - afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get()); + afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), bestAFguess); double reverseLog10PofNull = log10AlleleFrequencyPosteriors.get()[0]; - bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get()); double reverseLog10PofF = log10AlleleFrequencyPosteriors.get()[bestAFguess]; //System.out.println("reverseLog10PofNull=" + reverseLog10PofNull + ", reverseLog10PofF=" + reverseLog10PofF);