Have clip-vg trim flanks of big unaligned sequences - #1990
Merged
Conversation
Both aligners work by extending anchors, so an anchor in unique sequence beside
a satellite array gets extended into it, and inside the array the orthologous
alignment is undefined. Each haplotype takes some arbitrary but locally
high-scoring path for a while past the boundary, and they disagree with each
other. clip-vg's -u removes the runs longer than its threshold; what survives
is a tangle at every boundary, in runs too short for -u to reach.
clip-vg v1.2.2 adds -k/--flank for this, and this wires it up:
--clipFlank alongside --clip, how far to trim outward. This is the
control that matters -- inside a repeat the gate seldom
closes on its own, so the cap decides how much goes.
clipFlank config, default 50000, so it is on by default.
clipFlankThreshold config, negative to calibrate per graph, which is the
default and what you want.
One command-line option, the rest config, because only the cap is worth
reaching for from a command line. --clipFlank rejects a negative value and
rejects being combined with --clip 0; both would otherwise fall through the same
test as 0 and disable the trim silently, which next to clipFlankThreshold -- where
negative is the meaningful "calibrate" value -- is an easy thing to type.
The calibrated threshold decides how much sequence goes and clip-vg reports it
on stderr, which cactus_call only surfaces when a command fails, so clip_vg()
captures it into a flank-stats.tsv per chromosome. The threshold gets its own
column because that is the number worth reading down; the whole message sits
beside it so a rewording upstream costs the column rather than the record.
Validated whole-genome on HPRC v2.1 against CHM13, on the same chrom-alignments
as the release so only the clipping differs. Graph: nodes -1.53%, edges -1.98%,
length -1.49%, concentrated where it should be (chrY -9.4% of edges, the
acrocentrics next, chr6 -0.1%). VCF: raw -0.92%, final -0.61%. Accuracy
against dipcall truth over 231 HPRC samples, aardvark: F1 improves in 231 of
231, driven by recall.
Needs hal2vg v1.2.2, pinned here: v1.2.1 predates -k, so with clipFlank on by
default every run using downloaded binaries died at the clip phase on an unknown
option. All six hal2vg binaries move together, since they ship from one release.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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.
Right now, 10kb+ unaligned regions are clipped out by default in pangenome graphs (except in .full output). But sometimes, the flanks of these regions are poorly aligned (often extending into satellites). This PR adds logic to trim these flanks if they seem too sparsely aligned.
This is done with the new clip-vg -k option. From each edge of a clipped interval it walks outward node by node, up to clipFlank bp, scoring every unaligned base +(1−T) and every aligned base −T — so the running total over any stretch is u − L·T, positive exactly when that stretch's unaligned fraction exceeds T. It cuts at the furthest position where that running total was at its maximum, clipping everything between there and the original edge. Because it takes the maximum over the whole window rather than stopping at the first negative, it will cross an aligned stretch if enough unaligned sequence resumes beyond it.
This scheme is pretty conservative. Drops some censat sequence from hprc graphs, without much impact of vcf-related benchmarks or path/fragmentation stats. There may be room to be more aggressive later on.