Realignment filter #4698

davidbenjamin · 2018-04-24T16:20:17Z

Closes #2930.

@takutosato Here it is, finally. It's working as well as I can get it to work, which is decently but not a slam dunk.

droazen · 2018-04-24T19:10:52Z

src/main/java/org/broadinstitute/hellbender/engine/ReadsContext.java

+     * @param fragmentSize the maximum distance on either side of {@code interval} to look for mates.
+     * @return a map of reads ot their mates for all reads for which a mate could be found.
+     */
+    public Map<GATKRead, GATKRead> getReadToMateMap(final int fragmentSize) {


This might be more at home in ReadUtils rather than here, as it's a pretty specific method. You could add a generic query(interval) method to ReadsContext to enable this method to move to ReadUtils.

codecov-io · 2018-04-24T20:12:58Z

Codecov Report

Merging #4698 into master will increase coverage by 0.037%.
The diff coverage is 93.617%.

@@               Coverage Diff               @@
##              master     #4698       +/-   ##
===============================================
+ Coverage     80.111%   80.148%   +0.037%     
- Complexity     17442     17485       +43     
===============================================
  Files           1083      1083               
  Lines          63165     63313      +148     
  Branches       10193     10237       +44     
===============================================
+ Hits           50602     50744      +142     
+ Misses          8574      8573        -1     
- Partials        3989      3996        +7

Impacted Files	Coverage Δ	Complexity Δ
...ute/hellbender/utils/variant/GATKVCFConstants.java	`80% <ø> (ø)`	`4 <0> (ø)`	⬇️
...broadinstitute/hellbender/engine/ReadsContext.java	`100% <100%> (ø)`	`10 <1> (+1)`	⬆️
...a/org/broadinstitute/hellbender/utils/Trilean.java	`100% <100%> (ø)`	`3 <2> (+2)`	⬆️
...alignmentfilter/RealignmentArgumentCollection.java	`100% <100%> (ø)`	`1 <1> (?)`
...e/hellbender/utils/variant/GATKVCFHeaderLines.java	`99.301% <100%> (+0.005%)`	`10 <0> (ø)`	⬇️
...roadinstitute/hellbender/utils/read/ReadUtils.java	`79.81% <100%> (+0.592%)`	`196 <6> (+7)`	⬆️
...rs/realignmentfilter/FilterAlignmentArtifacts.java	`91.228% <91.228%> (ø)`	`19 <19> (?)`
...s/walkers/realignmentfilter/RealignmentEngine.java	`93.333% <93.333%> (ø)`	`35 <35> (?)`
...utils/smithwaterman/SmithWatermanIntelAligner.java	`50% <0%> (-30%)`	`1% <0%> (-2%)`
...adinstitute/hellbender/utils/spark/SparkUtils.java	`72.727% <0%> (-11.797%)`	`12% <0%> (-9%)`
... and 33 more

magicDGS

Just some random comments. I am trying to understand the algorithm to be able to apply the method too, so some of them might not make sense...

magicDGS · 2018-04-25T09:06:58Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+    public static final String MIN_MAPPING_QUALITY_LONG_NAME = "min-mapping-quality";
+    @Argument(fullName = MIN_MAPPING_QUALITY_LONG_NAME,
+            doc="Ignore reads with mapping quality less than this", optional=true)
+    private int minMappingQuality = DEFAULT_MIN_MAPPING_QUALITY;


This can be a default ReadFilter (MappingQualityReadFilter), which can be tweaked by the user as a plugin (delete the filter, change the minimum MQ, set a maximum one, etc). Otherwise, the filter can be applied twice!

Good point! Done.

magicDGS · 2018-04-25T09:08:56Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/realignmentfilter/Realigner.java

+        final List<BwaMemAlignment> alignments = aligner.alignSeqs(Arrays.asList(read), GATKRead::getBasesNoCopy).get(0);
+
+        final List<BwaMemAlignment> nonAltAlignments = alignments.size() == 1 ? alignments :
+                alignments.stream().filter(a -> a.getRefId() < numberOfRegularContigs).collect(Collectors.toList());


If there is only one alignment, it can still be alt, no?

Perhaps, although I haven't seen such a case yet. However, if there's only one alignment then we're not worried about a variant falsely appearing multi-mapped.

magicDGS · 2018-04-25T09:09:57Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/realignmentfilter/Realigner.java

+        } else if (alignments.get(0).getRefId() < 0) {
+            return new RealignmentResult(false, alignments);
+        } else if (alignments.size() == 1) {
+            return new RealignmentResult(true, alignments);


If there is only one alignment, is it really always in the supposed location?

Empirically, yes. I actually wrote a version of this tool that used a liftover chain file to map the hg38 alignment back to the original reference in order to check this, but it slowed the tool, complicated the code, made an additional resource file (which, unlike the hg38 index, depends on the original reference) necessary, and made no difference in the results.

Basically, hg19 is good enough that even if it misses some alignments its best alignment will always show up somewhere in the hg38 alignments.

This might not be the case for other reference genomes (e.g., non-model species). I thought that the filter can be used with other species too, which have much worse assemblies...

takutosato

Sorry this took a while. I get the sense that I still don't get the fundamental ideas of the realignment filter. In particular, the filter considers a realignment a success when the read realigns to a different chromosome - although I could be wrong.

I'd also be curious to hear more about the reasoning behind the filtering criteria, which I think goes as follows:

If more than 3 alt reads fail to realign, FILTER
Else, if we have 2 or more successful alignments, PASS
Else, FILTER if (# failed alignments > # success)

Could we discuss in person on Monday?

takutosato · 2018-04-25T18:29:54Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+ * <p>
+ *     The bam input to this tool should be the reassembly bamout produced by HaplotypeCaller or Mutect2 in the process of generating
+ *     the input callset.  The original bam will also work but might fail to filter some indels.  The reference passed with the -R argument
+ *     but be the reference to which the input bam was realigned.  This does not need to correspond to the reference of the BWAS-MEM


typo: but be -> must be?

typo: BWAS-MEM -> BWA-MEM

takutosato · 2018-04-26T15:12:51Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/realignmentfilter/Realigner.java

+        return checkAlignments(nonAltAlignments);
+    }
+
+    public RealignmentResult checkAlignments(final List<BwaMemAlignment> alignments) {


can be private

takutosato · 2018-04-26T15:16:41Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/realignmentfilter/Realigner.java

+        aligner.setFlagOption(BwaMemAligner.MEM_F_ALL);
+    }
+
+    public RealignmentResult realign(final GATKRead read) {


As suggested by IntelliJ, this method can be package-private (i.e. we can remove the public keyword). I don't know how big of a deal this is. One the one hand, it's good style to keep the visibility as tight as possible. On the other hand I get the sense that GATK code as a whole doesn't think it's that big of a deal

In the GATK I think we informally make anything available to a CLI public and use package-private largely for visible-for-testing methods.

I see. Thanks!

takutosato · 2018-04-26T15:32:22Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2IntegrationTest.java

@@ -68,6 +67,7 @@ public void testDreamTumorNormal(final File tumorBam, final String tumorSample,
        Utils.resetRandomGenerator();
        final File unfilteredVcf = createTempFile("unfiltered", ".vcf");
        final File filteredVcf = createTempFile("filtered", ".vcf");
+        final File alignmentFilteredVcf = createTempFile("alignment-filtered", ".vcf");


This file is not used

takutosato · 2018-04-26T15:40:43Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+    public void apply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext refContext, final FeatureContext fc) {
+        Trilean passesFilter = Trilean.UNKNOWN;
+
+        if (vc.getNAlleles() == 1) {


Does this expression ever evaluate to true?

Regardless, the ternary operator here will save a couple lines of code

Does this expression ever evaluate to true?

I think probably not for vcf in-spec, but I'm being defensive against non-spec vcfs and spec changes.

ternary. . .

done

takutosato · 2018-04-27T17:44:45Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+        }
+    }
+
+    private boolean supportsVariant(final GATKRead read, final VariantContext vc) {


It seems risky to me to have no unit tests for these helper methods - it's difficult to know whether they're working as expected.

You're right. I will be a good citizen and write some tests.

takutosato · 2018-04-27T17:47:26Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+                                .filter(m -> m.getRefId() == r.getRefId())
+                                .filter(m -> Math.abs(m.getRefStart() - r.getRefStart()) < realignmentArgumentCollection.maxReasonableFragmentLength)
+                                .filter(m -> SAMFlag.READ_REVERSE_STRAND.isSet(m.getSamFlag()) == SAMFlag.READ_REVERSE_STRAND.isSet(r.getSamFlag()))
+                                .map(m -> new Pair<BwaMemAlignment, BwaMemAlignment>(r,m)))


.map(m -> new Pair<BwaMemAlignment, BwaMemAlignment>(r,m))) -> .map(m -> new Pair<>(r,m)))

takutosato · 2018-04-27T19:07:16Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+                continue;
+            }
+
+            // First check the read itself.  If it maps uniquely we don't need to bother with the mate


If it maps uniquely we don't need to bother with the mate: is no longer accurate?

Good point. Deleted.

takutosato · 2018-04-27T19:41:50Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/realignmentfilter/Realigner.java

+        // TODO: we need to check that contig is the same or equivalent up to hg38 alt contig
+        // TODO: do this by looking at index.getReferenceContigNames() and bamHeader.getSequenceDictionary().getSequences()
+        // TODO: in IDE and seeing what the correspondence could be
+        // TODO: put in check that start position is within eg 10 Mb of original mapping


If the alt reads maps not exactly at the original location but within 10 Mb, we consider it a success because it could be the same location when lifted over to grch38?

That's the idea. This check is an edge case that might never actually arise because if the mapping is unique it's almost certainly in the right place. That's why I feel fine leaving this for later -- or never. But yes, the 10 MB is a liftover fudge factor. In practice I think the largest discrepancy between hg19 and hg38 coordinates is about 4 MB.

That might not be true for other assemblies (e.g., non-model species).

takutosato · 2018-04-27T19:45:37Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+                                .map(m -> new Pair<BwaMemAlignment, BwaMemAlignment>(r,m)))
+                        .collect(Collectors.toList());
+
+                if (plausiblePairs.size() <= 1) {


It seems we consider it successful when a read and its mate originally map to somewhere in chr3 but are both realigned uniquely to chr17

We do, but as I discussed above with @magicDGS this simply never happens in practice. Maybe if someone is using a really ancient reference.

Or any other species assembly. I would like to see this working for non-human data.

davidbenjamin · 2018-05-17T20:25:55Z

Back to @takutosato.

takutosato · 2018-05-18T17:49:55Z

src/main/java/org/broadinstitute/hellbender/utils/read/ReadUtils.java

+     * The algorithm is:
+     * 1) make two maps of read name --> read for reads overlapping {@code interval}, one for first-of-pair reads and one
+     *    for second-of-pair reads.
+     * 2) For all reads in an expanded interval padded by {@code fragmentSize} on both sides look for a read of the same


"read of the same that..." -> "read of the same name that..."

takutosato · 2018-05-18T20:39:17Z

.../org/broadinstitute/hellbender/tools/walkers/realignmentfilter/FilterAlignmentArtifacts.java

+                                .map(m -> new Pair<BwaMemAlignment, BwaMemAlignment>(r,m)))
+                        .collect(Collectors.toList());
+
+                if (plausiblePairs.size() <= 1) {


takutosato · 2018-05-18T21:21:35Z

...org/broadinstitute/hellbender/tools/walkers/realignmentfilter/RealignmentEngineUnitTest.java

+
+    @Test
+    public void testSupportsVariant() {
+        final GATKRead refRead = makeRead("refRead", 1, "ACGTACGT", "8M");


refRead is not used, but I think we should keep the ref sequence (extended to include "GGGGG") somehow, either as a comment or as an unused variable

I had meant to make sure that the ref read doesn't support any of the variants and forgot to do so. Fixed.

takutosato · 2018-05-18T21:28:18Z

...org/broadinstitute/hellbender/tools/walkers/realignmentfilter/RealignmentEngineUnitTest.java

+import java.util.Collection;
+import java.util.List;
+
+public class RealignmentEngineUnitTest {


Great tests!

takutosato · 2018-05-18T21:34:57Z

@davidbenjamin I think there was a typo in a comment but other than that looks great! Thanks for answering my questions - it didn't let me comment back (the reply boxes disappeared for some reason) but I went through them all.

Please merge when ready

davidbenjamin · 2018-05-22T13:42:19Z

@magicDGS you make some good points about other species. We can't devote any time to that right now, but I do think it's worth coming back to.

magicDGS · 2018-05-23T06:35:44Z

@davidbenjamin - can you please open an issue for my concerns and what is necessary for solving the issues? I might have a look at some point!

davidbenjamin · 2018-05-23T14:48:24Z

@magicDGS Done -- see issue #4803.

magicDGS · 2018-05-24T08:31:31Z

Thanks!

davidbenjamin added CLI Mutect labels Apr 24, 2018

davidbenjamin added this to the Mutect2 milestone Apr 24, 2018

davidbenjamin assigned takutosato Apr 24, 2018

davidbenjamin requested a review from takutosato April 24, 2018 16:20

davidbenjamin force-pushed the db_realignment_filter branch 3 times, most recently from c309c84 to 3eeef83 Compare April 24, 2018 19:05

droazen reviewed Apr 24, 2018

View reviewed changes

magicDGS reviewed Apr 25, 2018

View reviewed changes

takutosato requested changes Apr 27, 2018

View reviewed changes

takutosato assigned davidbenjamin and unassigned takutosato Apr 27, 2018

davidbenjamin force-pushed the db_realignment_filter branch 3 times, most recently from e6e5f7d to 9d603c7 Compare May 17, 2018 19:25

davidbenjamin assigned takutosato and unassigned davidbenjamin May 17, 2018

takutosato approved these changes May 18, 2018

View reviewed changes

Realignment filter

31c4f06

davidbenjamin force-pushed the db_realignment_filter branch from 9d603c7 to 31c4f06 Compare May 22, 2018 02:11

davidbenjamin merged commit ce0722a into master May 22, 2018

davidbenjamin deleted the db_realignment_filter branch May 22, 2018 13:39

Realignment filter #4698

Realignment filter #4698

Conversation

davidbenjamin commented Apr 24, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov-io commented Apr 24, 2018 • edited Loading

Codecov Report

magicDGS 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

takutosato 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

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davidbenjamin commented May 17, 2018

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

takutosato commented May 18, 2018

davidbenjamin commented May 22, 2018

magicDGS commented May 23, 2018

davidbenjamin commented May 23, 2018

magicDGS commented May 24, 2018

codecov-io commented Apr 24, 2018 •

edited

Loading