diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilterMutectCalls.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilterMutectCalls.java index 6c6bae4871b..a59512e5948 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilterMutectCalls.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilterMutectCalls.java @@ -104,7 +104,7 @@ public final class FilterMutectCalls extends MultiplePassVariantWalker { private static final int NUMBER_OF_LEARNING_PASSES = 2; @Override - protected int numberOfPasses() { return NUMBER_OF_LEARNING_PASSES + 1; } + protected int numberOfPasses() { return NUMBER_OF_LEARNING_PASSES + 2; } // {@code NUMBER_OF_LEARNING_PASSES} passes for learning, one for the threshold, and one for calling @Override public boolean requiresReference() { return true;} @@ -142,9 +142,9 @@ protected void nthPassApply(final VariantContext variant, final FeatureContext featureContext, final int n) { ParamUtils.isPositiveOrZero(n, "Passes must start at the 0th pass."); - if (n < NUMBER_OF_LEARNING_PASSES) { + if (n <= NUMBER_OF_LEARNING_PASSES) { filteringEngine.accumulateData(variant, referenceContext); - } else if (n == NUMBER_OF_LEARNING_PASSES) { + } else if (n == NUMBER_OF_LEARNING_PASSES + 1) { vcfWriter.add(filteringEngine.applyFiltersAndAccumulateOutputStats(variant, referenceContext)); } else { throw new GATKException.ShouldNeverReachHereException("This walker should never reach (zero-indexed) pass " + n); @@ -156,6 +156,10 @@ protected void afterNthPass(final int n) { if (n < NUMBER_OF_LEARNING_PASSES) { filteringEngine.learnParameters(); } else if (n == NUMBER_OF_LEARNING_PASSES) { + // it's important for filter parameters to stay the same and only learn the threshold in the final pass so that the + // final threshold used corresponds exactly to the filters + filteringEngine.learnThreshold(); + }else if (n == NUMBER_OF_LEARNING_PASSES + 1) { final Path filteringStats = IOUtils.getPath( filteringStatsOutput != null ? filteringStatsOutput : outputVcf + FILTERING_STATS_EXTENSION); diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java index c01324039e4..3a977b981b4 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java @@ -63,7 +63,7 @@ public Mutect2FilteringEngine(M2FiltersArgumentCollection MTFAC, final VCFHeader /** * Maximum probability that a potential variant is not a true somatic mutation. Variants with error probabilities - * below this threshold are called; variants with error probabilities above are filtered. + * at or below this threshold are called; variants with error probabilities above are filtered. */ public double getThreshold() { return thresholdCalculator.getThreshold(); } @@ -159,6 +159,11 @@ public void learnParameters() { filteringOutputStats.clear(); } + public void learnThreshold() { + thresholdCalculator.relearnThresholdAndClearAcumulatedProbabilities(); + filteringOutputStats.clear(); + } + /** * Create a filtered variant and record statistics for the final pass of {@link FilterMutectCalls} */ @@ -177,8 +182,9 @@ public VariantContext applyFiltersAndAccumulateOutputStats(final VariantContext } }); - // TODO: clarify this logic - if (errorProbability > EPSILON && errorProbability > getThreshold() - EPSILON) { + // error probability must exceed threshold, and just in case threshold is bad, probabilities close to 1 must be filtered + // and probabilities close to 0 must not be filtered + if (( errorProbability > Math.min(1 - EPSILON, Math.max(EPSILON, getThreshold())))) { vcb.filter(entry.getKey().filterName()); } }