Hard filters in M2 for mitochondria #5842

meganshand · 2019-03-28T15:03:49Z

This adds a hard filter for low variant allele fraction calls. We will not turn this on by default in any of our pipelines, but it will give users an easy option to filter everything below a certain VAF that they don't care about.

It also adds a hard filter for low alt depth calls based on a threshold from the median autosomal coverage (that must be supplied as an argument). It takes the cutoff from a Poisson with a mean of 1.5 * median coverage (to account for NuMTs with 3 copies in the autosome) and is tuned to catch 99% of the false positives (which we know will also catch lots of true positives).

It also removes the Polymorphic NUMT annotation (since that's basically what's going into the filter at this point).

codecov-io · 2019-03-28T15:07:24Z

Codecov Report

Merging #5842 into master will increase coverage by 0.004%.
The diff coverage is 94.444%.

@@              Coverage Diff               @@
##              master    #5842       +/-   ##
==============================================
+ Coverage     87.056%   87.06%   +0.004%     
- Complexity     32157    32178       +21     
==============================================
  Files           1978     1979        +1     
  Lines         147433   147498       +65     
  Branches       16206    16215        +9     
==============================================
+ Hits          128350   128412       +62     
- Misses         13171    13172        +1     
- Partials        5912     5914        +2

Impacted Files	Coverage Δ	Complexity Δ
...itute/hellbender/tools/walkers/mutect/Mutect2.java	`80.952% <ø> (-0.866%)`	`22 <0> (-1)`
...der/tools/walkers/mutect/M2ArgumentCollection.java	`90% <ø> (ø)`	`14 <0> (ø)`	⬇️
...ute/hellbender/utils/variant/GATKVCFConstants.java	`75% <ø> (ø)`	`4 <0> (ø)`	⬇️
...r/tools/walkers/mutect/Mutect2IntegrationTest.java	`87.776% <ø> (-0.021%)`	`91 <0> (ø)`
...lkers/mutect/filtering/Mutect2FilteringEngine.java	`97.115% <100%> (+0.057%)`	`43 <0> (ø)`	⬇️
.../mutect/filtering/M2FiltersArgumentCollection.java	`93.75% <100%> (+0.417%)`	`6 <0> (ø)`	⬇️
...kers/mutect/filtering/MinAlleleFractionFilter.java	`100% <100%> (ø)`	`7 <7> (?)`
...e/hellbender/utils/variant/GATKVCFHeaderLines.java	`94.886% <100%> (+0.029%)`	`11 <0> (ø)`	⬇️
...alkers/mutect/filtering/PolymorphicNuMTFilter.java	`88.235% <88.235%> (ø)`	`9 <9> (?)`
...lbender/utils/variant/GATKVariantContextUtils.java	`87.309% <0%> (-0.306%)`	`244% <0%> (-2%)`
... and 10 more

droazen · 2019-03-28T15:31:56Z

@meganshand Is this something you are targeting for the current release? Our current plan is to release towards the end of the day today.

meganshand · 2019-03-28T15:36:18Z

@droazen while that would be amazing, I don't expect @ldgauthier or @davidbenjamin to drop everything to review this (unless one of you happens to have some free time today?). Don't let this hold back a release.

davidbenjamin · 2019-03-28T15:40:06Z

I've got some time to try.

davidbenjamin

Done with review, back to @meganshand.

davidbenjamin · 2019-03-28T15:45:36Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/MinAfFilter.java

+
+    @Override
+    public boolean isArtifact(final VariantContext vc, final Mutect2FilteringEngine filteringEngine) {
+        MutableBoolean b = new MutableBoolean(false);


This could be cast as an anyMatch:

return vc.getGenotypes().stream().filter(filteringEngine::isTumor) .filter(g -> g.hasExtendedAttribute(GATKVCFConstants.ALLELE_FRACTION_KEY)) .anyMatch(g -> { final double[] alleleFractions = GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(g, GATKVCFConstants.ALLELE_FRACTION_KEY, () -> null, 1.0); double max = alleleFractions[0]; for (int i=0; i<alleleFractions.length; i++) { //Only include allele fractions of non symbolic alleles for GVCFs if (!(vc.hasSymbolicAlleles() && i == alleleFractions.length - 1)) { max = alleleFractions[i] > max ? alleleFractions[i] : max; } } return max < minAf; });

davidbenjamin · 2019-03-28T15:49:37Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/MinAfFilter.java

+                .filter(g -> g.hasExtendedAttribute(GATKVCFConstants.ALLELE_FRACTION_KEY))
+                .forEach(g -> {
+                    final double[] alleleFractions = GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(g, GATKVCFConstants.ALLELE_FRACTION_KEY, () -> null, 1.0);
+                    double max = alleleFractions[0];


The maximization can be shortened, and max should be final:

final int numRealAlleles = vc.hasSymbolicAlleles() ? alleleFractions.length - 1 : alleleFractions.length; final double max = IntStream.range(0, numRealAlleles).mapToDouble(a -> alleleFractions[a]).max();

davidbenjamin · 2019-03-28T15:50:21Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/MinAfFilter.java

+
+import java.util.*;
+
+public class MinAfFilter extends HardFilter {


Expand name to MinAlleleFractionFilter

davidbenjamin · 2019-03-28T15:51:17Z

...java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/PolymorphicNuMTFilter.java

+public class PolymorphicNuMTFilter extends HardFilter {
+    private static final double LOWER_BOUND_PROB = .01;
+    private static final double MULTIPLE_COPIES_MULTIPLIER = 1.5;
+    private int maxAltDepthCutoff;


can this be final?

davidbenjamin · 2019-03-28T15:52:35Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/MinAfFilter.java

+    }
+
+    @Override
+    public ErrorType errorType() { return ErrorType.ARTIFACT; }


It's your judgment call, but I'm torn between ErrorType.ARTIFACT and ErrorType.SEQUENCING.

Since these are mostly mapping artifacts (from NuMTs) I was leaning towards ARTIFACT, but maybe let me know if you still think this should be SEQUENCING.

Oh I see your comment below, should these be NON_SOMATIC? Really low VAF things includes sequencing errors, mapping errors, etc.

davidbenjamin · 2019-03-28T15:53:12Z

...java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/PolymorphicNuMTFilter.java

+    }
+
+    @Override
+    public ErrorType errorType() { return ErrorType.ARTIFACT; }


This should be ErrorType.NON_SOMATIC because it's a mapping error of real DNA.

davidbenjamin · 2019-03-28T15:54:00Z

...java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/PolymorphicNuMTFilter.java

+
+    @Override
+    public boolean isArtifact(final VariantContext vc, final Mutect2FilteringEngine filteringEngine) {
+        MutableBoolean b = new MutableBoolean(false);


Similar comments as above above using anyMatch on the Stream and making the max code more concise.

davidbenjamin · 2019-03-28T15:54:21Z

...java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/PolymorphicNuMTFilter.java

+    public boolean isArtifact(final VariantContext vc, final Mutect2FilteringEngine filteringEngine) {
+        MutableBoolean b = new MutableBoolean(false);
+        vc.getGenotypes().stream().filter(filteringEngine::isTumor)
+                .filter(g -> g.hasAD())


Use a method reference Genotype::hasAD

davidbenjamin · 2019-03-28T15:55:36Z

src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java

@@ -173,6 +172,8 @@ their names (or descriptions) depend on some threshold.  Those filters are not i
    public final static String STRICT_STRAND_BIAS_FILTER_NAME =               "strict_strand";
    public final static String N_RATIO_FILTER_NAME =                           "n_ratio";
    public final static String CHIMERIC_ORIGINAL_ALIGNMENT_FILTER_NAME =       "numt_chimera"; //mitochondria
+    public final static String ALLELE_FRACTION_FILTER_NAME =                   "low_allele_fraction";


For beautiful formatting of FilterMutectCalls' filter analysis output, could you makes the filter names more than 8 and less than 16 characters?

meganshand · 2019-03-28T19:37:29Z

Thanks @davidbenjamin! This looks much cleaner now. Still holding out a tiny bit of hope for merging before the release, if you have a chance to look at the new changes.

davidbenjamin

👍 and good luck merging in time!

droazen · 2019-03-28T21:01:19Z

@meganshand Latest build failed with:

org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.testWithAllAnnotations FAILED
    java.lang.AssertionError: Iterators differ at element [9]: FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification"> != FORMAT=<ID=NUMT,Number=1,Type=String,Description="Potentially a polymorphic NuMT false positive rather than a real mitochondrial variant."> expected [FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">] but found [FORMAT=<ID=NUMT,Number=1,Type=String,Description="Potentially a polymorphic NuMT false positive rather than a real mitochondrial variant.">]
        at org.testng.Assert.fail(Assert.java:93)
        at org.testng.Assert.failNotEquals(Assert.java:512)
        at org.testng.Assert.assertEqualsImpl(Assert.java:134)
        at org.testng.Assert.assertEquals(Assert.java:610)
        at org.testng.Assert.assertEquals(Assert.java:578)
        at org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.assertHeadersMatch(VariantAnnotatorIntegrationTest.java:66)
        at org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.runVariantAnnotatorAndAssertSomething(VariantAnnotatorIntegrationTest.java:92)
        at org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.assertVariantContextsMatch(VariantAnnotatorIntegrationTest.java:55)
        at org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.assertVariantContextsMatch(VariantAnnotatorIntegrationTest.java:49)
        at org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorIntegrationTest.testWithAllAnnotations(VariantAnnotatorIntegrationTest.java:224)
Exception in thread "Thread-5" java.nio.file.ClosedFileSystemException
	at com.google.common.jimfs.FileSystemState.checkOpen(FileSystemState.java:64)
	at com.google.common.jimfs.JimfsFileStore.lookUp(JimfsFileStore.java:141)
	at com.google.common.jimfs.FileSystemView.lookUp(FileSystemView.java:123)
	at com.google.common.jimfs.FileSystemView.lookUpWithLock(FileSystemView.java:112)
	at com.google.common.jimfs.FileSystemView.readAttributes(FileSystemView.java:771)
	at com.google.common.jimfs.JimfsFileSystemProvider.readAttributes(JimfsFileSystemProvider.java:346)
	at java.nio.file.Files.readAttributes(Files.java:1737)
	at java.nio.file.FileTreeWalker.getAttributes(FileTreeWalker.java:219)
	at java.nio.file.FileTreeWalker.visit(FileTreeWalker.java:276)
	at java.nio.file.FileTreeWalker.walk(FileTreeWalker.java:322)
	at java.nio.file.Files.walkFileTree(Files.java:2662)
	at java.nio.file.Files.walkFileTree(Files.java:2742)
	at htsjdk.samtools.util.IOUtil.recursiveDelete(IOUtil.java:1344)
	at org.broadinstitute.hellbender.utils.io.IOUtils.deleteRecursively(IOUtils.java:1061)
	at org.broadinstitute.hellbender.utils.io.DeleteRecursivelyOnExitPathHook.runHooks(DeleteRecursivelyOnExitPathHook.java:56)
	at java.lang.Thread.run(Thread.java:748)
Results: FAILURE (1604 tests, 1603 successes, 1 failures, 0 skipped)

meganshand requested review from davidbenjamin and ldgauthier March 28, 2019 15:04

droazen assigned davidbenjamin Mar 28, 2019

davidbenjamin requested changes Mar 28, 2019

View reviewed changes

davidbenjamin approved these changes Mar 28, 2019

View reviewed changes

Adding hard filter for low VAF and polymorphic NuMTs

2738424

meganshand force-pushed the ms_new_filters branch from 2efdf6a to 2738424 Compare March 28, 2019 19:54

Removing header line of deleted annotation

d39fb0b

meganshand merged commit ea3032d into master Mar 28, 2019

meganshand deleted the ms_new_filters branch March 28, 2019 22:33

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Hard filters in M2 for mitochondria #5842

Hard filters in M2 for mitochondria #5842

meganshand commented Mar 28, 2019

codecov-io commented Mar 28, 2019 •

edited

Loading

droazen commented Mar 28, 2019

meganshand commented Mar 28, 2019

davidbenjamin commented Mar 28, 2019

davidbenjamin left a comment

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

meganshand Mar 28, 2019

meganshand Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

davidbenjamin Mar 28, 2019

meganshand commented Mar 28, 2019

davidbenjamin left a comment

droazen commented Mar 28, 2019


		import java.util.*;

		public class MinAfFilter extends HardFilter {

Hard filters in M2 for mitochondria #5842

Hard filters in M2 for mitochondria #5842

Conversation

meganshand commented Mar 28, 2019

codecov-io commented Mar 28, 2019 • edited Loading

Codecov Report

droazen commented Mar 28, 2019

meganshand commented Mar 28, 2019

davidbenjamin commented Mar 28, 2019

davidbenjamin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

meganshand commented Mar 28, 2019

davidbenjamin left a comment

Choose a reason for hiding this comment

droazen commented Mar 28, 2019

codecov-io commented Mar 28, 2019 •

edited

Loading