diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilterMutectCalls.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilterMutectCalls.java index d3194b4ed3e..04c6affbf29 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilterMutectCalls.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilterMutectCalls.java @@ -22,8 +22,6 @@ import java.io.File; import java.util.Optional; -import java.util.ArrayList; -import java.util.List; import java.util.Set; import java.util.stream.Collectors; @@ -84,9 +82,7 @@ public final class FilterMutectCalls extends TwoPassVariantWalker { private Mutect2FilteringEngine filteringEngine; - private List firstPassFilterResults; - - private Mutect2FilterSummary stats; + private FilteringFirstPass filteringFirstPass; @Override public void onTraversalStart() { @@ -109,7 +105,7 @@ public void onTraversalStart() { final Optional normalSample = normalSampleHeaderLine == null ? Optional.empty() : Optional.of(normalSampleHeaderLine.getValue()); filteringEngine = new Mutect2FilteringEngine(MTFAC, tumorSample, normalSample); - firstPassFilterResults = new ArrayList<>(); + filteringFirstPass = new FilteringFirstPass(tumorSample); } @Override @@ -120,18 +116,18 @@ public Object onTraversalSuccess() { @Override public void firstPassApply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext refContext, final FeatureContext fc) { final FilterResult filterResult = filteringEngine.calculateFilters(MTFAC, vc, Optional.empty()); - firstPassFilterResults.add(filterResult); + filteringFirstPass.add(filterResult, vc); } @Override protected void afterFirstPass() { - stats = filteringEngine.calculateFilterStats(firstPassFilterResults, MTFAC.maxFalsePositiveRate); - Mutect2FilterSummary.writeM2FilterSummary(stats, MTFAC.mutect2FilteringStatsTable); + filteringFirstPass.learnModelForSecondPass(MTFAC.maxFalsePositiveRate); + filteringFirstPass.writeM2FilterSummary(MTFAC.mutect2FilteringStatsTable); } @Override public void secondPassApply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext refContext, final FeatureContext fc) { - final FilterResult filterResult = filteringEngine.calculateFilters(MTFAC, vc, Optional.of(stats)); + final FilterResult filterResult = filteringEngine.calculateFilters(MTFAC, vc, Optional.of(filteringFirstPass)); final VariantContextBuilder vcb = new VariantContextBuilder(vc); vcb.filters(filterResult.getFilters()); diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPass.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPass.java new file mode 100644 index 00000000000..89498fd1bf4 --- /dev/null +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPass.java @@ -0,0 +1,213 @@ +package org.broadinstitute.hellbender.tools.walkers.mutect; + + +import htsjdk.variant.variantcontext.Genotype; +import htsjdk.variant.variantcontext.VariantContext; +import org.apache.commons.lang3.tuple.ImmutablePair; +import org.apache.commons.lang3.tuple.Pair; +import org.broadinstitute.hellbender.exceptions.UserException; +import org.broadinstitute.hellbender.utils.Utils; +import org.broadinstitute.hellbender.utils.param.ParamUtils; +import org.broadinstitute.hellbender.utils.tsv.DataLine; +import org.broadinstitute.hellbender.utils.tsv.TableColumnCollection; +import org.broadinstitute.hellbender.utils.tsv.TableWriter; +import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants; + +import java.io.File; +import java.io.IOException; +import java.util.*; + +/** + * Stores the results of the first pass of {@link FilterMutectCalls}, a purely online step in which each variant is + * not "aware" of other variants, and learns various global properties necessary for a more refined second step. + */ +public class FilteringFirstPass { + final List filterResults; + final Map> filteredPhasedCalls; + final Map filterStats; + final String tumorSample; + boolean readyForSecondPass; + + public FilteringFirstPass(final String tumorSample) { + filterResults = new ArrayList<>(); + filteredPhasedCalls = new HashMap<>(); + filterStats = new HashMap<>(); + readyForSecondPass = false; + this.tumorSample = tumorSample; + } + + public boolean isReadyForSecondPass() { return readyForSecondPass; } + + public FilterStats getFilterStats(final String filterName){ + Utils.validateArg(filterStats.containsKey(filterName), "invalid filter name: " + filterName); + return filterStats.get(filterName); + } + + public boolean isOnFilteredHaplotype(final VariantContext vc, final int maxDistance) { + + final Genotype tumorGenotype = vc.getGenotype(tumorSample); + + if (!hasPhaseInfo(tumorGenotype)) { + return false; + } + + final String pgt = (String) tumorGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY, ""); + final String pid = (String) tumorGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY, ""); + final int position = vc.getStart(); + + final Pair filteredCall = filteredPhasedCalls.get(pid); + if (filteredCall == null) { + return false; + } + + // Check that vc occurs on the filtered haplotype + return filteredCall.getLeft().equals(pgt) && Math.abs(filteredCall.getRight() - position) <= maxDistance; + } + + public void add(final FilterResult filterResult, final VariantContext vc) { + filterResults.add(filterResult); + final Genotype tumorGenotype = vc.getGenotype(tumorSample); + + if (!filterResult.getFilters().isEmpty() && hasPhaseInfo(tumorGenotype)) { + final String pgt = (String) tumorGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY, ""); + final String pid = (String) tumorGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY, ""); + final int position = vc.getStart(); + filteredPhasedCalls.put(pid, new ImmutablePair<>(pgt, position)); + } + } + + public void learnModelForSecondPass(final double requestedFPR){ + final double[] readOrientationPosteriors = getFilterResults().stream() + .filter(r -> r.getFilters().isEmpty()) + .mapToDouble(r -> r.getReadOrientationPosterior()) + .toArray(); + + final FilterStats readOrientationFilterStats = calculateThresholdForReadOrientationFilter(readOrientationPosteriors, requestedFPR); + filterStats.put(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, readOrientationFilterStats); + readyForSecondPass = true; + } + + /** + * + * Compute the filtering threshold that ensures that the false positive rate among the resulting pass variants + * will not exceed the requested false positive rate + * + * @param posteriors A list of posterior probabilities, which gets sorted + * @param requestedFPR We set the filtering threshold such that the FPR doesn't exceed this value + * @return + */ + public static FilterStats calculateThresholdForReadOrientationFilter(final double[] posteriors, final double requestedFPR){ + ParamUtils.isPositiveOrZero(requestedFPR, "requested FPR must be non-negative"); + final double thresholdForFilteringNone = 1.0; + final double thresholdForFilteringAll = 0.0; + + Arrays.sort(posteriors); + + final int numPassingVariants = posteriors.length; + double cumulativeExpectedFPs = 0.0; + + for (int i = 0; i < numPassingVariants; i++){ + final double posterior = posteriors[i]; + + // One can show that the cumulative error rate is monotonically increasing in i + final double expectedFPR = (cumulativeExpectedFPs + posterior) / (i + 1); + if (expectedFPR > requestedFPR){ + return i > 0 ? + new FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, posteriors[i-1], + cumulativeExpectedFPs, i-1, cumulativeExpectedFPs/i, requestedFPR) : + new FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, thresholdForFilteringAll, + 0.0, 0, 0.0, requestedFPR); + } + + cumulativeExpectedFPs += posterior; + } + + // If the expected FP rate never exceeded the max tolerable value, then we can let everything pass + return new FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, thresholdForFilteringNone, + cumulativeExpectedFPs, numPassingVariants, cumulativeExpectedFPs/numPassingVariants, requestedFPR); + } + + public static boolean hasPhaseInfo(final Genotype genotype) { + return genotype.hasExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY) && genotype.hasExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY); + } + + public List getFilterResults() { + return filterResults; + } + + public static class FilterStats { + private final String filterName; + private final double threshold; + private final double expectedNumFPs; + private final int numPassingVariants; + private final double expectedFPR; + private final double requestedFPR; + + public FilterStats(final String filterName, final double threshold, final double expectedNumFPs, + final int numPassingVariants, final double expectedFPR, final double requestedFPR){ + this.filterName = filterName; + this.threshold = threshold; + this.expectedNumFPs = expectedNumFPs; + this.numPassingVariants = numPassingVariants; + this.expectedFPR = expectedFPR; + this.requestedFPR = requestedFPR; + } + + public String getFilterName() { return filterName; } + + public double getExpectedNumFPs() { return expectedNumFPs; } + + public int getNumPassingVariants() { return numPassingVariants; } + + public double getThreshold() { return threshold; } + + public double getExpectedFPR() { return expectedFPR; } + + public double getRequestedFPR() { return requestedFPR; } + + } + + private enum M2FilterStatsTableColumn { + FILTER_NAME("filter_name"), + THRESHOLD("threshold"), + EXPECTED_FALSE_POSITIVES("expected_fps"), + EXPECTED_FALSE_POSITIVE_RATE("expected_fpr"), + REQUESTED_FALSE_POSITIVE_RATE("requested_fpr"), + NUM_PASSING_VARIANTS("num_passing_variants"); + + private String columnName; + + M2FilterStatsTableColumn(final String columnName) { + this.columnName = columnName; + } + + @Override + public String toString() { return columnName; } + + public static final TableColumnCollection COLUMNS = new TableColumnCollection((Object[]) values()); + } + + private static class Mutect2FilterStatsWriter extends TableWriter { + private Mutect2FilterStatsWriter(final File output) throws IOException { + super(output, M2FilterStatsTableColumn.COLUMNS); + } + + @Override + protected void composeLine(final FilterStats stats, final DataLine dataLine) { + dataLine.set(M2FilterStatsTableColumn.FILTER_NAME.toString(), stats.getFilterName()) + .set(M2FilterStatsTableColumn.THRESHOLD.toString(), stats.getThreshold()) + .set(M2FilterStatsTableColumn.EXPECTED_FALSE_POSITIVES.toString(), stats.getExpectedNumFPs()) + .set(M2FilterStatsTableColumn.EXPECTED_FALSE_POSITIVE_RATE.toString(), stats.getExpectedFPR()) + .set(M2FilterStatsTableColumn.REQUESTED_FALSE_POSITIVE_RATE.toString(), stats.getRequestedFPR()) + .set(M2FilterStatsTableColumn.NUM_PASSING_VARIANTS.toString(), stats.getNumPassingVariants()); + } + } + + public void writeM2FilterSummary(final File outputTable) { + try (Mutect2FilterStatsWriter writer = new Mutect2FilterStatsWriter(outputTable)) { + writer.writeAllRecords(filterStats.values()); + } catch (IOException e) { + throw new UserException(String.format("Encountered an IO exception while writing to %s.", outputTable), e); + } + } +} diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2FiltersArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2FiltersArgumentCollection.java index 59a764a1758..dfa83e36a19 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2FiltersArgumentCollection.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2FiltersArgumentCollection.java @@ -28,6 +28,7 @@ public class M2FiltersArgumentCollection extends AssemblyBasedCallerArgumentColl public static final String UNIQUE_ALT_READ_COUNT_LONG_NAME = "unique-alt-read-count"; public static final String TUMOR_SEGMENTATION_LONG_NAME = "tumor-segmentation"; public static final String ORIENTATION_BIAS_FDR_LONG_NAME = "orientation-bias-fdr"; // FDR = false discovery rate + public static final String MAX_DISTANCE_TO_FILTERED_CALL_ON_SAME_HAPLOTYPE_LONG_NAME = "distance-on-haplotype"; public static final String FILTERING_STATS_LONG_NAME = "stats"; @@ -124,4 +125,8 @@ public class M2FiltersArgumentCollection extends AssemblyBasedCallerArgumentColl public File mutect2FilteringStatsTable = new File("Mutect2FilteringStats.tsv"); + @Argument(fullName = MAX_DISTANCE_TO_FILTERED_CALL_ON_SAME_HAPLOTYPE_LONG_NAME, optional = true, doc = "On second filtering pass, variants with same PGT and PID tags as a filtered variant within this distance are filtered.") + public int maxDistanceToFilteredCallOnSameHaplotype = 100; + + } diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilterSummary.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilterSummary.java deleted file mode 100644 index d2dc30ada24..00000000000 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilterSummary.java +++ /dev/null @@ -1,129 +0,0 @@ -package org.broadinstitute.hellbender.tools.walkers.mutect; - - -import org.broadinstitute.hellbender.exceptions.UserException; -import org.broadinstitute.hellbender.utils.Utils; -import org.broadinstitute.hellbender.utils.tsv.DataLine; -import org.broadinstitute.hellbender.utils.tsv.TableColumnCollection; -import org.broadinstitute.hellbender.utils.tsv.TableWriter; - -import java.io.File; -import java.io.IOException; -import java.util.Collection; -import java.util.HashMap; -import java.util.Map; -import java.util.Optional; - -/** - * - */ -public class Mutect2FilterSummary { - final Map map; - - public Mutect2FilterSummary(){ - map = new HashMap<>(); - } - - public void addNewFilterStats(final String filterName, final FilterStats filterStats){ - map.put(filterName, filterStats); - } - - public FilterStats getFilterStats(final String filterName){ - Utils.validateArg(map.containsKey(filterName), "invalid filter name: " + filterName); - return map.get(filterName); - } - - public Collection getAllFilterStats(){ - return map.values(); - } - - public static class FilterStats { - private final String filterName; - private final double threshold; - private final double expectedNumFPs; - private final int numPassingVariants; - private final double expectedFPR; - private final double requestedFPR; - - public FilterStats(final String filterName, final double threshold, final double expectedNumFPs, - final int numPassingVariants, final double expectedFPR, final double requestedFPR){ - this.filterName = filterName; - this.threshold = threshold; - this.expectedNumFPs = expectedNumFPs; - this.numPassingVariants = numPassingVariants; - this.expectedFPR = expectedFPR; - this.requestedFPR = requestedFPR; - } - - public String getFilterName() { - return filterName; - } - - public double getExpectedNumFPs() { - return expectedNumFPs; - } - - public int getNumPassingVariants() { - return numPassingVariants; - } - - public double getThreshold() { - return threshold; - } - - public double getExpectedFPR() { - return expectedFPR; - } - - public double getRequestedFPR() { - return requestedFPR; - } - - } - - private static class Mutect2FilterStatsWriter extends TableWriter { - private Mutect2FilterStatsWriter(final File output) throws IOException { - super(output, M2FilterStatsTableColumn.COLUMNS); - } - - @Override - protected void composeLine(final FilterStats stats, final DataLine dataLine) { - dataLine.set(M2FilterStatsTableColumn.FILTER_NAME.toString(), stats.getFilterName()) - .set(M2FilterStatsTableColumn.THRESHOLD.toString(), stats.getThreshold()) - .set(M2FilterStatsTableColumn.EXPECTED_FALSE_POSITIVES.toString(), stats.getExpectedNumFPs()) - .set(M2FilterStatsTableColumn.EXPECTED_FALSE_POSITIVE_RATE.toString(), stats.getExpectedFPR()) - .set(M2FilterStatsTableColumn.REQUESTED_FALSE_POSITIVE_RATE.toString(), stats.getRequestedFPR()) - .set(M2FilterStatsTableColumn.NUM_PASSING_VARIANTS.toString(), stats.getNumPassingVariants()); - } - } - - public static void writeM2FilterSummary(final Mutect2FilterSummary summary, final File outputTable) { - try (Mutect2FilterStatsWriter writer = new Mutect2FilterStatsWriter(outputTable)) { - writer.writeAllRecords(summary.getAllFilterStats()); - } catch (IOException e) { - throw new UserException(String.format("Encountered an IO exception while writing to %s.", outputTable), e); - } - } - - private enum M2FilterStatsTableColumn { - FILTER_NAME("filter_name"), - THRESHOLD("threshold"), - EXPECTED_FALSE_POSITIVES("expected_fps"), - EXPECTED_FALSE_POSITIVE_RATE("expected_fpr"), - REQUESTED_FALSE_POSITIVE_RATE("requested_fpr"), - NUM_PASSING_VARIANTS("num_passing_variants"); - - private String columnName; - - M2FilterStatsTableColumn(final String columnName) { - this.columnName = columnName; - } - - @Override - public String toString() { - return columnName; - } - - public static final TableColumnCollection COLUMNS = new TableColumnCollection((Object[]) values()); - } -} diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngine.java index 8c0f7303307..8cbc884d070 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngine.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngine.java @@ -8,10 +8,7 @@ import org.broadinstitute.hellbender.tools.walkers.annotator.*; import org.broadinstitute.hellbender.tools.walkers.contamination.ContaminationRecord; import org.broadinstitute.hellbender.tools.walkers.contamination.MinorAlleleFractionRecord; -import org.broadinstitute.hellbender.utils.GATKProtectedVariantContextUtils; -import org.broadinstitute.hellbender.utils.IndexRange; -import org.broadinstitute.hellbender.utils.MathUtils; -import org.broadinstitute.hellbender.utils.QualityUtils; +import org.broadinstitute.hellbender.utils.*; import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants; import java.util.*; import java.util.stream.Collectors; @@ -295,7 +292,7 @@ private void applyDuplicatedAltReadFilter(final M2FiltersArgumentCollection MTFA } } - private void applyReadOrientationFilter(final VariantContext vc, final FilterResult filterResult, final Optional filterSummary){ + private void applyReadOrientationFilter(final VariantContext vc, final FilterResult filterResult, final Optional firstPass){ if (! vc.isSNP()){ return; } @@ -308,12 +305,12 @@ private void applyReadOrientationFilter(final VariantContext vc, final FilterRes final double artifactPosterior = GATKProtectedVariantContextUtils.getAttributeAsDouble(tumorGenotype, GATKVCFConstants.ROF_POSTERIOR_KEY, -1.0); - if (! filterSummary.isPresent()) { + if (! firstPass.isPresent()) { // During first pass we simply collect the posterior artifact probabilities filterResult.setReadOrientationPosterior(artifactPosterior); return; } else { - final double threshold = filterSummary.get().getFilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME).getThreshold(); + final double threshold = firstPass.get().getFilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME).getThreshold(); if (artifactPosterior > threshold){ filterResult.addFilter(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME); @@ -321,9 +318,17 @@ private void applyReadOrientationFilter(final VariantContext vc, final FilterRes } } + private void applyFilteredHaplotypeFilter(final M2FiltersArgumentCollection MTFAC, final VariantContext vc, final FilterResult filterResult, final Optional firstPass){ + if ( firstPass.isPresent() && firstPass.get().isOnFilteredHaplotype(vc, MTFAC.maxDistanceToFilteredCallOnSameHaplotype)){ + filterResult.addFilter(GATKVCFConstants.BAD_HAPLOTYPE_FILTER_NAME); + } + } + public FilterResult calculateFilters(final M2FiltersArgumentCollection MTFAC, final VariantContext vc, - final Optional filterStats) { + final Optional firstPass) { + firstPass.ifPresent(ffp -> Utils.validate(ffp.isReadyForSecondPass(), "First pass information has not been processed into a model for the second pass.")); final FilterResult filterResult = new FilterResult(); + applyFilteredHaplotypeFilter(MTFAC, vc, filterResult, firstPass); applyInsufficientEvidenceFilter(MTFAC, vc, filterResult); applyClusteredEventFilter(vc, filterResult); applyDuplicatedAltReadFilter(MTFAC, vc, filterResult); @@ -340,7 +345,7 @@ public FilterResult calculateFilters(final M2FiltersArgumentCollection MTFAC, fi applyReadPositionFilter(MTFAC, vc, filterResult); // The following filters use the information gathered during the first pass - applyReadOrientationFilter(vc, filterResult, filterStats); + applyReadOrientationFilter(vc, filterResult, firstPass); return filterResult; } @@ -348,56 +353,4 @@ private int[] getIntArrayTumorField(final VariantContext vc, final String key) { return GATKProtectedVariantContextUtils.getAttributeAsIntArray(vc.getGenotype(tumorSample), key, () -> null, 0); } - /** - * - * Compute the filtering threshold that ensures that the false positive rate among the resulting pass variants - * will not exceed the requested false positive rate - * - * @param posteriors A list of posterior probabilities, which gets sorted - * @param requestedFPR We set the filtering threshold such that the FPR doesn't exceed this value - * @return - */ - public Mutect2FilterSummary.FilterStats calculateThresholdForReadOrientationFilter(final double[] posteriors, final double requestedFPR){ - final double thresholdForFilteringNone = 1.0; - final double thresholdForFilteringAll = 0.0; - - Arrays.sort(posteriors); - - final int numPassingVariants = posteriors.length; - double cumulativeExpectedFPs = 0.0; - - for (int i = 0; i < numPassingVariants; i++){ - final double posterior = posteriors[i]; - - // One can show that the cumulative error rate is monotonically increasing in i - final double expectedFPR = (cumulativeExpectedFPs + posterior) / (i + 1); - if (expectedFPR > requestedFPR){ - return i > 0 ? - new Mutect2FilterSummary.FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, posteriors[i-1], - cumulativeExpectedFPs, i-1, cumulativeExpectedFPs/i, requestedFPR) : - new Mutect2FilterSummary.FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, thresholdForFilteringAll, - 0.0, 0, 0.0, requestedFPR); - } - - cumulativeExpectedFPs += posterior; - } - - // If the expected FP rate never exceeded the max tolerable value, then we can let everything pass - return new Mutect2FilterSummary.FilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, thresholdForFilteringNone, - cumulativeExpectedFPs, numPassingVariants, cumulativeExpectedFPs/numPassingVariants, requestedFPR); - } - - public Mutect2FilterSummary calculateFilterStats(final List filterResults, final double requestedFPR){ - final Mutect2FilterSummary filterSummary = new Mutect2FilterSummary(); - - final double[] readOrientationPosteriors = filterResults.stream() - .filter(r -> r.getFilters().isEmpty()) - .mapToDouble(r -> r.getReadOrientationPosterior()) - .toArray(); - - final Mutect2FilterSummary.FilterStats readOrientationFilterStats = calculateThresholdForReadOrientationFilter(readOrientationPosteriors, requestedFPR); - filterSummary.addNewFilterStats(GATKVCFConstants.READ_ORIENTATION_ARTIFACT_FILTER_NAME, readOrientationFilterStats); - - return filterSummary; - } } diff --git a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java index ad7195c3d5c..ef27535522a 100644 --- a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java +++ b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java @@ -158,6 +158,7 @@ their names (or descriptions) depend on some threshold. Those filters are not i public final static String READ_POSITION_FILTER_NAME = "read_position"; public final static String CONTAMINATION_FILTER_NAME = "contamination"; public final static String READ_ORIENTATION_ARTIFACT_FILTER_NAME = "read_orientation_artifact"; + public final static String BAD_HAPLOTYPE_FILTER_NAME = "bad_haplotype"; public static final List MUTECT_FILTER_NAMES = Arrays.asList(STR_CONTRACTION_FILTER_NAME, PON_FILTER_NAME, CLUSTERED_EVENTS_FILTER_NAME, TUMOR_LOD_FILTER_NAME, GERMLINE_RISK_FILTER_NAME, @@ -165,7 +166,7 @@ their names (or descriptions) depend on some threshold. Those filters are not i MEDIAN_BASE_QUALITY_FILTER_NAME, MEDIAN_MAPPING_QUALITY_FILTER_NAME, MEDIAN_FRAGMENT_LENGTH_DIFFERENCE_FILTER_NAME, READ_POSITION_FILTER_NAME, CONTAMINATION_FILTER_NAME, DUPLICATED_EVIDENCE_FILTER_NAME, - READ_ORIENTATION_ARTIFACT_FILTER_NAME); + READ_ORIENTATION_ARTIFACT_FILTER_NAME, BAD_HAPLOTYPE_FILTER_NAME); // Symbolic alleles public final static String SYMBOLIC_ALLELE_DEFINITION_HEADER_TAG = "ALT"; diff --git a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java index b6690d06101..504a61c44a0 100644 --- a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java +++ b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java @@ -60,6 +60,7 @@ private static void addFilterLine(final VCFFilterHeaderLine line) { addFilterLine(new VCFFilterHeaderLine(CONTAMINATION_FILTER_NAME, "contamination")); addFilterLine(new VCFFilterHeaderLine(DUPLICATED_EVIDENCE_FILTER_NAME, "evidence for alt allele is overrepresented by apparent duplicates")); addFilterLine(new VCFFilterHeaderLine(READ_ORIENTATION_ARTIFACT_FILTER_NAME, "orientation bias detected by the orientation bias mixture model")); + addFilterLine(new VCFFilterHeaderLine(BAD_HAPLOTYPE_FILTER_NAME, "Variant near filtered variant on same haplotype.")); addFormatLine(new VCFFormatHeaderLine(ALLELE_BALANCE_KEY, 1, VCFHeaderLineType.Float, "Allele balance for each het genotype")); diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngineUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPassUnitTest.java similarity index 80% rename from src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngineUnitTest.java rename to src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPassUnitTest.java index 93286017461..b688a988a85 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2FilteringEngineUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/FilteringFirstPassUnitTest.java @@ -8,14 +8,7 @@ import java.util.Optional; -public class Mutect2FilteringEngineUnitTest extends BaseTest { - Mutect2FilteringEngine engine; - - @BeforeSuite - public void init(){ - engine = new Mutect2FilteringEngine(new M2FiltersArgumentCollection(), "badTumor", Optional.empty()); - } - +public class FilteringFirstPassUnitTest extends BaseTest { @DataProvider(name = "falsePositiveRateData") public Object[][] makeFalsePositiveRateData() { return new Object[][]{ @@ -33,7 +26,7 @@ public Object[][] makeFalsePositiveRateData() { public void testCalculateThresholdForReadOrientationFilter(final double[] posteriors, final double maxErrorRate, final double expectedThreshold){ - final Mutect2FilterSummary.FilterStats stats = engine.calculateThresholdForReadOrientationFilter(posteriors, maxErrorRate); + final FilteringFirstPass.FilterStats stats = FilteringFirstPass.calculateThresholdForReadOrientationFilter(posteriors, maxErrorRate); Assert.assertEquals(stats.getThreshold(), expectedThreshold); }