Trimming support in gaf overlap filter - #2008
Merged
Merged
Conversation
…records GAFOverlapFilterTrim, off by default, passes gaffilter -t so a record that loses a query overlap gives up only the contested span rather than the whole record. Whole-record deletion also destroys the uncontested flanks, which on HPRC v2.1 chr15 at 464 haplotypes is about half of the 1.373 Gb it removes. The trim takes its own GAFOverlapFilterTrimMinLengthRatio, defaulting to 0, rather than inheriting GAFOverlapFilterMinLengthRatio. That guard exists because deleting a record is expensive and a trim is not, so with trimming on it only leaves small conflicts between long contigs unadjudicated: measured on a 12-sample chr15 graph, keeping it at 0.25 while trimming doubly places 409,584 bp on one segmental-duplication pair that 0 does not. Flipping the one attribute named Trim therefore selects the configuration that was measured, not a worse neighbour of it. GAFOverlapFilterTrimRescueWeak (gaffilter -R) is off and should stay off absent evidence this run cannot provide. A cut can leave an unaligned hole with alignment on both sides, and closing it means handing the span to one claimant -- but no such choice can meet the bar GAFOverlapFilterRatio sets, because a hole only forms where every record spanning it lost to a record it does not dominate. On the same graph, declining costs 7.97 Mb of recovered path but avoids 25,299 nodes and 121,960 bp of new graph sequence: the recovery that needs no decision lands at 1,023 bp per new node and makes the graph 1,267 bp shorter, while the guessed part lands at 315 bp per new node and accounts for essentially all of the growth. Measured end to end, 12-sample chr15, against the current default: +10,103,364 bp of haplotype path, +20 path fragments, graph 1,267 bp shorter, no haplotype losing sequence and the reference untouched. The -l file gaffilter needs is written by gaf2unstable in full before it emits its first GAF line, and gaffilter reads all of its input before opening it, so there is no race in the pipe -- but only in that order, which is why -l is passed the same path gaf2unstable -o writes. NOTE: build-tools/downloadPangenomeTools still pins a cactus-gfa-tools commit without -t. With the attribute off nothing changes, and turning it on against that pin fails loudly rather than silently, but the pin has to be bumped to the gaffilter commit once it is pushed, in a follow-up, before this is usable. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MCSd5hqBNbBk11pcbC53Wx
Passes gaffilter -e. The bases abutting an overlap are the least certain part of the alignment, so with trimming on the cut extends a little past the contested span on each side. 5000 is half the cactus-graphmap-join clip threshold, and that is the point: where only the losing record pulls back the gap it leaves is exactly this size, so cactus carries it instead of clipping it. At 10000 it would sit on the boundary. Cost measured on HPRC v2.1 at 464 haplotypes, as a fraction of the block length the trim keeps: chr15 0.8%, chr9 1.7%, chr21 5.2%, chr20 9.6% -- largest where the absolute amount recovered is smallest. Holes barely move (chr15 1132 -> 1138, chr9 1008 -> 1012). Still gated behind GAFOverlapFilterTrim, which remains off. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MCSd5hqBNbBk11pcbC53Wx
…e records Passes gaffilter -Q. A record that loses a query overlap and is itself poorly placed keeps the behaviour it has today, deleted whole, rather than contributing its uncontested flanks to the graph: losing an overlap was always evidence about the record and not only about the overlapping part. 20 sits above minMAPQ=5, so it is strictly more conservative than filter_paf, and well below the mapq-60 mass that carries the benefit. Block length kept, at 464 haplotypes: chr15 -0.04%, chr9 -1.68%, chr20 -0.97%, chr21 -19.17%. Still gated behind GAFOverlapFilterTrim, which remains off. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MCSd5hqBNbBk11pcbC53Wx
GAFOverlapFilterTrim needs a gaffilter built with -t, and the previous pin (b37287b1) predates it: turning the attribute on against that build failed every minigraph_map_one job with "invalid option -- 't'". The attribute is off by default, so nothing changed for anyone, but the config documented a capability the shipped binary did not have. a32588c is on cactus-gfa-tools overlap-work, not yet merged to main. It stays valid as long as that branch is not squashed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MCSd5hqBNbBk11pcbC53Wx
Of the six the trim originally added, three could not earn their place. GAFOverlapFilterTrimMinGap and GAFOverlapFilterTrimRescueWeak are a dead pair. --close-holes is off and cannot be justified -- a hole only forms where every record spanning it lost to one it does not dominate, so no claimant to it can meet the bar GAFOverlapFilterRatio sets -- and with it off, -g provably does not change the output: identical md5 over -g 0, 1000, 10000 and 1000000 on the 12-sample chr15 GAF. Both remain available as gaffilter flags for a hand-run. GAFOverlapFilterTrimMinLengthRatio was a knob with one defensible setting. The guard it mirrors exists because deleting a record is expensive, which a trim is not; trimming at 0.25 doubly places 409,584 bp on one segmental-duplication pair that 0 does not, and nothing favours a middle value. It is now just 0, with the reason in a comment, so the attribute named Trim still selects the configuration that was measured. Left: GAFOverlapFilterTrim, GAFOverlapFilterTrimEdge, GAFOverlapFilterTrimMinMAPQ, each with a measured effect. Output is unchanged -- the command this now builds is byte-identical in result to the one the measured arms ran. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MCSd5hqBNbBk11pcbC53Wx
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Note: this is off-by-default. But am merging it now as I've re-invented flavours of it several times while trying to clean up complex regions, and I think this is the most promising incarnation.
The root issue in many cases is poor inversion handling in the base sv graph (tackled in #2006), but this work should provide a complimentary angle. Needs full-genome testing before activiating by default.
Basically instead of removing "dominated" mappings in presence of query overlaps, it cuts ends (plus a trim factor).