Improved allele clustering - #4998
Merged
Merged
Conversation
Deconstructor::get_alleles spells a star allele as the literal "*" for an empty Traversal, then treats that marker as if it were sequence in three places: the "substitution" length fold, reverse_complement_in_place, and the prev_char padding. complement['*'] is 'N' (src/utility.cpp), so at a reverse-oriented site the spanning-deletion marker silently became a real ambiguous base. The four observed spellings were "*", "N", "AN" and "A*". "A*" is invalid VCF -- htslib only recognizes an overlapping-deletion ALT when the field is exactly "*". Bare "N" is worse because it is silently valid: vg construct warns and skips a variant whose ALT is "*", so a deconstruct -> construct round trip turned "skip this variant" into a real N node with nothing downstream objecting. Identify stars by traversal emptiness rather than by the string, and skip them in all three steps. Also guard the empty-reference case: the star's length-1 "*" used to be what forced padded (indel) form at a star-only insertion site, and without that the record would emit REF ".". Removing the star from the substitution fold makes the spelling of the real alleles independent of whether -R was passed; records containing a star can change position as a result, since add_variant sorts on (contig, position). Adds five fixtures under test/nesting/ covering forward, reverse, padded, reverse+padded, MNP and empty-REF sites, plus a sweep asserting that every allele whose AT entry is "." is spelled exactly "*". Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
-L/--cluster and --cluster-post were parsed, range-validated, cross-validated and passed into the FlowCaller constructor, but never read anywhere. The option had been a silent no-op since v1.73.0, and 18_vg_call.t asserted that "clustering does not change number of variant sites" -- a test that passed precisely because nothing happened. Implement -L as what it is actually wanted for: collapsing a 1/2 call whose two alleles are effectively the same into 1/1 with a single ALT. The parameters mirror vg deconstruct exactly -- -L/--cluster F plus --cluster-min-len N, same defaults, same site gate -- and the metric is inherited rather than reimplemented, by calling deconstruct's own cluster_traversals(). So the endpoint pruning, the weighted-Jaccard and the >= comparison are identical: on the same graph both tools flip at the same threshold (verified by bisection at 5/7 in the new tests). The merge runs at the very end of VCFOutputCaller::emit_variant, after update_vcf_info and after flatten_common_allele_ends. Merging earlier is unsound: AD/DP are recomputed inside update_vcf_info over site_traversals, so an absorbed allele's reads vanish rather than pooling, and PoissonSupportSnarlCaller::genotype_likelihood charges a Poisson penalty over the traversals *not* in the genotype, so removing one perturbs every genotype score. Running after flattening keeps the surviving allele's string, POS and REF byte-identical to a run without -L. The --cluster-min-len site gate reads the alleles the record emits, not the traversal finder's candidate list. The finder returns up to max_yens_traversals (50) speculative paths, most of which never become an allele, so gating on those would let a branch with no reads and no AT entry decide whether merging happens. Clustering is done in descending allele-depth order so that each cluster's head -- the surviving allele, and the one MAT's Jaccard is measured against -- is its best-supported member. The finder's own ranking switches to length-weighted average flow on snarls past the average-support threshold, and can put a short lightly-supported allele ahead of a long heavily-supported one; merging into that head would emit the minority sequence as a homozygous call carrying the pooled depth. AD is summed onto the survivor, GL is folded by max over the merged genotype classes, MAD is recomputed, and the merge is recorded in a new MAT info field so a merged 1/1 is distinguishable from a real hom-alt. GQ/GP/DP/QUAL are deliberately left alone: they describe the site and come from the caller's own CallInfo rather than the emitted GL. Merging is ALT-vs-ALT only. vg deconstruct folds near-reference alleles into the reference cluster and drops the record; doing that here would turn a het call into no call at all, so it is deliberately not copied. Two divergences from deconstruct, both intentional: out-of-range -L is an error rather than clamped (-L 5 is a plausible typo for -L 0.5), and -L is rejected with -v, -G/-T and -B, which never reach emit_variant or whose statistics a merge would invalidate. --cluster-post is removed: with merging always post-genotyping it has no meaning, and its documented "output grouping only" was never implementable. Also fixes two stale error messages that referred to the legacy caller as -L. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…nterior
--cluster-min-len promised "50 = SVs only" but measured the snarl traversal
interior, which is not the size of the variant. A 1bp SNP sitting inside a
121bp snarl passed the gate and was clustered away:
vg call g.vg -k g.pack -p x x 101 T C,G 1/2
vg call g.vg -k g.pack -p x -L 0.6 --cluster-min-len 50
x 101 T G 1/1
The gate bisected at 121/122 -- the interior -- nowhere near the 1bp record.
The failure was positively correlated with the gate: a large shared interior
is both what raises the Jaccard and what opens the gate, so --cluster-min-len
offered small variants essentially no protection. vg deconstruct had the same
defect; on that graph it deleted the SNP record outright at 50.
Both tools now gate on CORE LENGTH: the longest allele after stripping the
prefix and suffix common to every non-"*" allele. One helper,
VCFOutputCaller::allele_core_length, used by both.
Core length is the right measure precisely because it is invariant to how much
shared context a tool keeps in its strings. vg call flattens its alleles down
to an anchor base and vg deconstruct emits the whole snarl interior, so a raw
string length answers differently for the same variant; the core does not.
That invariance also makes the gate independent of -a, which flattens
differently again.
Three consequences, all wanted. The anchor base that flattening must leave on
every indel is a shared prefix and is stripped, so a 49bp indel measures 49 and
not 50. The reference allele participates, because a pure deletion emits
|ALT|=1 and a maximum over ALTs alone would never fire for a deletion of any
size. And "*" is excluded, which also neutralises flatten_common_allele_ends
being a no-op whenever a "*" is present.
In vg deconstruct the allele strings do not exist at gate time -- get_alleles
runs after clustering -- so the gate rebuilds the allele set the record would
emit un-clustered. That is exact, not approximate: the only transforms
get_alleles applies afterwards are reverse-complementing every allele and
prepending a common anchor base, and core length is invariant to both. The old
interior-length loop survives as an exact pre-filter, since core length can
never exceed the longest interior, so nothing gets slower.
The change is one-directional: core <= max interior, so sites stop clustering
but none start. yeast100 output is byte-identical with the gate active.
Note core length measures the variant's span, not the size of any single event
inside it: a haplotype differing from the reference at two bases 59bp apart has
a core length of 60. The alternative -- the difference in allele lengths --
measures 0 for a 60bp substitution, which is worse.
Also corrects three stale comments left by the preceding commit: the
set_allele_merge doc and the deconstruct cross-reference, both falsified by
this change, and two misattributions (GL is written by
PoissonSupportSnarlCaller::update_vcf_info, not genotype; -B does reach
emit_variant and is rejected for a different reason than -v and -G/-T).
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Traversal clustering could not merge a pure deletion with anything. In
cluster_traversals each traversal became a multiset of its interior handles,
except that a 2-visit traversal -- a pure deletion, straight from snarl start
to snarl end -- kept BOTH BOUNDARY handles instead:
int64_t first = trav.size() == 2 ? 0 : 1;
int64_t last = trav.size() == 2 ? trav.size() : trav.size() - 1;
Every longer traversal kept only its interior, so the two sets were disjoint by
construction, the weighted Jaccard was structurally 0, and no threshold could
merge them. A 60bp deletion and a 59bp deletion of the same site failed to
merge at every -L down to 0.0.
That made -L's advertised use case inert: deletions are the dominant SV
representation in a vg graph, and --cluster-min-len is documented as the SV
knob. The standing "// todo: does jaccard properly handle empty sets?" beside
that code was pointing at the same thing.
A pure deletion now prunes to the empty set, which is what it is -- it carries
no sequence -- and the pair is scored against the site instead:
sim(deletion, X) = 1 - |X| / max(|site|, |X|)
the fraction of the site both alleles delete. 59 of 60bp gives 59/60.
The site is used ONLY when one allele is a pure deletion. Scaling every pair
by it would make two unrelated alleles look similar merely because they sit in
a big snarl; test/nesting/inv60_in_2kb.gfa and unrelated10_in_2kb.gfa pin that,
and both fail under a site-scaled-everywhere metric. A consequence worth
stating: every pair of traversals with three or more visits each scores exactly
as before, so no existing threshold changes meaning.
sim(A,A) is 1 for every input including two identical pure deletions, where the
arithmetic is 0/0. Returning NaN there would fail both the > and the >=
comparisons in cluster_traversals and make the traversal permanently
unclusterable -- the same class of bug being fixed.
cluster_traversals takes the site reference as a new trailing parameter with no
default, so a caller that omits it fails to compile rather than silently
reverting to the old behaviour. All three callers supply it: vg deconstruct
its reference traversal, vg call the VCF reference (nullptr for a
NestedFlowCaller child-snarl Visit, which has no handle representation and
falls back to pairwise scoring), vg simplify its chosen reference candidate.
vg simplify is the only consumer that mutates the graph, so a merged deletion
now removes the absorbed allele's nodes and edges: 4 nodes/5 edges to 3/3 on
the new fixture. It is gated behind -L < 1.0, which is off by default, and
that path had no -L coverage at all before this commit.
Adds a unit test on the similarity function itself. This defect survived
several reviews because the metric was only ever exercised through two layers
of caller, where a similarity of 0 is indistinguishable from "these alleles
really are unrelated".
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Fixes from a review of the preceding commits. No change to the merge or clustering logic itself. Restore the --bottom-up + -L rejection. It was removed on the argument that it had been written for a design that was never implemented, without checking whether -L actually works in that mode. It does not: NestedFlowCaller represents a child snarl as a Visit carrying a Snarl rather than a node, those have no handle for the clusterer, and the merge therefore does nothing at exactly the nested sites --bottom-up exists for -- while the MAT header advertises the feature. An accepted-but-inert option is the defect this feature was written to remove. Move the -L validations above the graph load. Every condition is pure command-line state, but they ran after loading and overlaying the graph, so a typo'd "-L 5" cost 59.6s and 4.0GB on chr22 before being rejected. Now 0.00s. Reject -L nan. It passed the range check, because both nan < 0.0 and nan > 1.0 are false, and then silently disabled merging -- precisely what the range check exists to prevent. -L inf was already rejected. Warn that -L has no effect at ploidy 1, where a genotype can never have the two distinct called ALTs a merge needs. Accept --cluster-post again, warning that it is deprecated and ignored. It shipped as a documented no-op through v1.76, so removing it outright turned any pipeline carrying it into an unrecognized-option failure. Remove after a release. Say "length-weighted Jaccard" in both tools' helptext. The metric weights by node length: alleles sharing half their nodes can score 0.968, so a user tuning F against the unweighted meaning lands far off. Use one parser for the AD field. It was being read three ways in one function, including vg::parse<double>, which exits rather than returning -- inside an OpenMP region, where that abandons whatever other threads have buffered. Comment corrections, all verified against the code: the parse helper's rationale named the wrong alternative (a non-exiting 2-argument vg::parse does exist; the real reason is that it can throw); "nothing gets slower" was false, the gate builds one interior string per traversal against one per cluster in get_alleles, roughly 2x on a many-haplotype site; the GL layout is written by PoissonSupportSnarlCaller::update_vcf_info, not genotype; the flatten no-op has a different mechanism with -a than without; "exact pre-filter" is one-sided, so "conservative"; and a third stale message describing --legacy as -L. Test fixes. The POS assertion compared a column the awk filter had already pinned, so it asserted 9 == 9. The GL assertion counted fields, which a spec-ordered fold would also satisfy; it now pins the values. The MAT header was only asserted absent without -L, never present with it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
None of -G, --bottom-up or -R/--ploidy-regex is exercised anywhere in the test
suite, which is how all three of these survived. All are pre-existing; none is
introduced by the -L work on this branch.
-G segfaulted on any graph. GAFOutputCaller::emit_gaf_variant builds gt_travs,
a NEW vector holding one entry per called allele, then passed the original
ref_trav_idx -- an index into the full traversal list -- straight through to
emit_gaf_traversals, which used it to index the small vector. It also pushed
travs[allele] for every genotype entry, including the negative star and missing
markers, which address nothing. Remap the reference index into gt_travs and
skip the markers; bounds-check the reference index in emit_gaf_traversals too,
since three other call sites pass it directly.
--bottom-up aborted on any nested graph whose reference traversal crosses a
child snarl. NestedFlowCaller plants Visits carrying a Snarl rather than a
node, and PoissonSupportSnarlCaller::update_vcf_info called
graph.get_handle(node_id) on them with node_id 0. Skip those visits; they
contribute no length of their own.
-R/--ploidy-regex accepted any ploidy, bypassing the {1,2} check -d gets, and
reached the caller as an unsupported ploidy. vg call -R 'x:3' crashed where
vg call -d 3 correctly errors. Validate the rule at parse time.
Adds the first tests for all three modes. Each fails against the unfixed code:
reverting the three changes reddens exactly tests 67, 68, 83, 85 and 86.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…n modes Follow-up to a review of this branch. The MAT info field described its value as OLD>NEW:JACCARD, but the number is not a Jaccard coefficient when a pure deletion is involved -- that is the whole point of the preceding commit, which renamed the function accordingly. On the new deletion fixture MAT reports 0.983 where the Jaccard is 0. The field is now OLD>NEW:SIMILARITY. The helptexts keep "length-weighted Jaccard", which is accurate for every pair that carries sequence; vg simplify's helptext said "(handle) Jaccard coefficient" and now matches the other two. Reject --bottom-up with -T and -G. NestedFlowCaller represents a child snarl as a Visit carrying a Snarl rather than a node; -T feeds those to to_mapping, which asserts, and -G emits a GAF header with no records. Both are pre-existing. Rejecting matches how the same Visits are handled for -L rather than leaving an accepted-but-inert option, which is the defect this branch exists to remove. Reject --legacy with -L. LegacyCaller has its own traversal finder and support model and has never been exercised with the merge; it was the one caller missing from the rejection list. The ploidy-1 warning tested -d only, so a haploid ploidy set through -R/--ploidy-regex slipped past it silently. Test corrections. del60_vs_snp.gfa is a 60bp substitution, not the "1bp SNP" two descriptions claimed; adds del1_vs_snp1.gfa, which covers the small-variant case none of the fixtures reached -- a 1bp deletion and a 1bp SNP score 0 and never merge, since deleting a base and substituting it are different events. The vg simplify assertions named a mechanism that does not occur: the merge removes an edge, and the node count falls because unchop then fuses the chain. The edge assertion is the one that pins the merge. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Correctness:
- vg simplify chose ref_trav_candidates[0] -- the alphabetically-first spanning
path -- where it meant the -P reference. At a top-level snarl every spanning
path is first seen there and so ties with the reference at rank 0, leaving the
name comparison to decide. A -P path now wins that tie outright. Two
consequences, one of which has nothing to do with -L:
- That traversal seeds the set of nodes and edges vg simplify KEEPS, so
whenever a haplotype sorted ahead of the reference, the reference's own
nodes were deleted and its path left in fragments. On test/yeast.vg,
`vg simplify -a small -P S288C` with no -L at all lost 4.5% of the
reference; it now preserves all 16 contigs and 12,157,149 bp
byte-for-byte, while still taking the graph from 546,797 to 354,513 nodes.
- It is also the scale a pure deletion is measured against, which made the
new -L merging depend on path names: del59_vs_del60.gfa and
simplify_del_absorbs.gfa are structurally identical and differ only in
their path names, and only the latter merged. The site is taken from a
path that was in the -P set, with nullptr (pairwise scoring, which merges
less) when none spans the snarl.
- Reject -L with star alleles in BOTH tools: vg call -Y and vg deconstruct -R.
A "*" in a child record means "an upstream deletion covers this site", and
clustering can absorb the allele that deletion came from -- which one survives
is decided by traversal order, so by path names -- leaving the "*" referring to
nothing in the file. That is a malformed record rather than a lossy one, which
is what separates it from ordinary nested clustering: there a clustered parent
deliberately disagrees with its own child records, giving the collapsed view of
a large variant while the children keep the precise one. That is allowed in
both tools and now tested, including the gref workflow, where the interior of
an insertion the reference bypasses becomes its own contig. vg deconstruct
-a -R -L emitted the malformed form before this branch.
- Move the -R/--ploidy-regex {1,2} check to where the rule is applied. Checking
at parse time rejected rules matching none of the called contigs, so
-R 'chrY:0' -p x went from exit 0 on v1.76.1 to exit 1.
- Hoist the site length out of weighted_traversal_similarity. It was re-summed
on every pairwise comparison, on the default vg deconstruct path: 2.32s ->
0.32s on a 10k-node reference interior with pure deletions present, matching
the no-deletion control.
- --cluster-min-len now defaults to 50 rather than 0 in both tools. The
similarity runs over the whole snarl interior, so at 0 it is dominated by
shared flanks: on core_snp_in_flanks.gfa a het T>C / T>G scores 0.984 and
vg deconstruct -L 0.9 folds both alleles into the reference cluster and drops
the record entirely. The gate only applies with -L, which is off by default,
so no default run changes. The "no effect without -L" warning now keys on
whether the flag was passed, or it would fire on every invocation.
Documentation:
- All three -L helptexts said "length-weighted Jaccard"; for a pure deletion the
compared value is 1 - |X|/max(|site|,|X|), which is the case -L exists for.
- Drop the stale one-liner left stacked above weighted_traversal_similarity's
doc block, and the stale "jaccard" in deconstructor.{hpp,cpp} and
traversal_clusters.hpp.
- graph_caller.cpp claimed AD is Number=R (it is Number=.) and that the merge
keeps DP == sum(AD) exact (not a vg invariant; the true claim is that sum(AD)
is unchanged), and that both tools gate the same variant the same way (the
gated sets differ). deconstructor.cpp's cost estimate was inverted.
- Trim the -L validation block and the unit-test file's framing.
Tests: 18_vg_call.t 158 -> 167, 26_deconstruct.t 109 -> 116, 43_vg_simplify.t
11 -> 16. Merge-mechanics assertions on small fixtures now pass
--cluster-min-len 0 explicitly, since their subject is the metric rather than
the gate; the new default is pinned separately in both tools. dladder had
become vacuous for del1_vs_snp1 under the new default (the gate decided a 1bp
site before the metric was consulted) and now measures the metric alone. New
simplify cases pin the site selection with a fixture whose alt paths sort before
the reference, and pin that the -P reference keeps its full length under the -m
length filter when a haplotype sorts ahead of it. A new gref case pins that -L merges a 121bp insertion at the
default gate while the nested record on the gref fragment stays untouched. Both "-L 1.0 is byte-identical to no -L" assertions were
tautologies -- there is no "was -L given" flag -- and now also assert the
absence of the MAT header and of TS/TL.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
No textual conflicts. One adaptation was needed, folded in here so the merge commit itself is green: 54bfd0f ("Spell deconstruct's snarl IDs in reference orientation, as call already does") flips v.id for a site the reference traverses backwards, so the child snarl in nesting/nested_snp_in_del_rev.gfa and nesting/nested_del_star_indel_rev.gfa is now spelled <5<2 rather than >2>5. Two star-allele assertions look records up by that ID and so found nothing. Only the key moved -- the records still read ID=<5<2 POS=3 REF=A ALT=T,* ID=<5<2 POS=2 REF=CAAAA ALT=CT,* i.e. exactly the "*" the star fix on this branch exists to produce, neither reverse-complemented into N nor padded. The two lookups are keyed on <5<2 and a comment records why the reverse fixtures differ from their forward twins.
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
2 participants
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.
Changelog Entry
To be copied to the draft changelog by merger:
vg deconstruct -Loption ported tovg callto cluster big alleles. Both now support deletions. Remains experimental.Description