From 2cc95e07a50c409fba620236e7d357a0f3debb25 Mon Sep 17 00:00:00 2001 From: David Benjamin Date: Mon, 19 Aug 2019 15:39:40 -0400 Subject: [PATCH 1/2] fixed Mutect2 bug that overfiltered by one variant --- .../walkers/mutect/filtering/FilterMutectCalls.java | 10 +++++++--- .../mutect/filtering/Mutect2FilteringEngine.java | 12 +++++++++--- 2 files changed, 16 insertions(+), 6 deletions(-) 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..9ab66219c1c 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; } // {@coode 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()); } } From c690880f126811a491242c7354e38d4b37cfe7a5 Mon Sep 17 00:00:00 2001 From: David Benjamin Date: Tue, 20 Aug 2019 22:50:01 -0400 Subject: [PATCH 2/2] whoops --- .../tools/walkers/mutect/filtering/FilterMutectCalls.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 9ab66219c1c..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 + 2; } // {@coode NUMBER_OF_LEARNING_PASSES} passes for learning, one for the threshold, and one for calling + 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;}