diff --git a/.gitignore b/.gitignore index 463350a3..683f43d7 100644 --- a/.gitignore +++ b/.gitignore @@ -23,4 +23,7 @@ Thumbs.db __pycache__/ ehthumbs.db -/debug \ No newline at end of file +/debug + +# compiled typst notes (source .typ is tracked) +/notes/*.pdf diff --git a/notes/assets/n00/dual_benchmark.csv b/notes/assets/n00/dual_benchmark.csv new file mode 100644 index 00000000..9cc35639 --- /dev/null +++ b/notes/assets/n00/dual_benchmark.csv @@ -0,0 +1,21 @@ +dataset,variant,seconds,n_paths,n_blocks,pangenome_bp,core_bp,core_blocks +example.fa,nodual,0.0,2,1,1000,1000,1 +example.fa,dual,0.0,2,1,1000,1000,1 +russian_doll_plasmids.fa.gz,nodual,0.2,4,11,50790,33610,5 +russian_doll_plasmids.fa.gz,dual,0.2,4,11,50790,33610,5 +flu-h3.fa,nodual,0.2,51,1,1737,1737,1 +flu-h3.fa,dual,0.2,51,1,1737,1737,1 +flu-h1.fa,nodual,0.8,171,10,2332,762,1 +flu-h1.fa,dual,0.9,171,10,2332,762,1 +ges-1.fa,nodual,7.8,33,218,596005,1016,1 +ges-1.fa,dual,8.2,33,218,596005,1016,1 +mpox.fa,nodual,19.3,13,63,207072,165129,19 +mpox.fa,dual,21.9,13,63,207072,165129,19 +sc2.fa,nodual,14.4,169,404,145717,1094,3 +sc2.fa,dual,17.0,169,404,145717,1094,3 +campylobacter-3.fa,nodual,17.3,3,512,3445626,102659,49 +campylobacter-3.fa,dual,21.8,3,514,3445993,102867,50 +klebs.fa.gz,nodual,142,9,1381,7644985,4458249,268 +klebs.fa.gz,dual,157,9,1372,7644453,4457866,268 +ecoli.fa.gz,nodual,157,10,2914,7830290,3782006,498 +ecoli.fa.gz,dual,187,10,2908,7827250,3782120,498 diff --git a/notes/assets/n00/dual_quality.png b/notes/assets/n00/dual_quality.png new file mode 100644 index 00000000..901a53b5 Binary files /dev/null and b/notes/assets/n00/dual_quality.png differ diff --git a/notes/assets/n00/dual_runtime.png b/notes/assets/n00/dual_runtime.png new file mode 100644 index 00000000..1857f028 Binary files /dev/null and b/notes/assets/n00/dual_runtime.png differ diff --git a/notes/n00_build_order_dependence.typ b/notes/n00_build_order_dependence.typ new file mode 100644 index 00000000..c4f48215 --- /dev/null +++ b/notes/n00_build_order_dependence.typ @@ -0,0 +1,773 @@ +#set text(font: "New Computer Modern", size: 11pt) +#set page(margin: 2cm) +#set par(justify: true) +#set heading(numbering: "1.1") +#show link: set text(fill: blue.darken(20%)) +#show raw.where(block: true): set block( + fill: luma(245), + inset: 8pt, + radius: 3pt, + width: 100%, +) + += Order dependence in `pangraph build` + +#emph[Branch `debug/order-dependence`. Investigation of the report in `tmp/pangraph_order_bug/`.] + +== Summary + +`pangraph build` produces a materially different graph depending only on the order the input +FASTA files are listed on the command line --- 452 to 509 blocks for the same three genomes, +a #sym.tilde 12% spread. The reported hypothesis was that graph merging follows argument order +rather than the guide tree. That is #emph[not] what happens: the guide tree is honoured exactly. + +The real cause is that block identifiers are assigned from the input position, those identifiers +are handed to `minimap2` as sequence names, and `minimap2` is invoked with `-X`, which keeps only +one direction of each pairwise alignment --- chosen by #emph[string comparison of those names]. +Argument order therefore decides, for every pair of blocks, which one is the query and which one +is the reference. That choice is not neutral: the alignment is asymmetric, and downstream the +reference is privileged when choosing the consensus that the merged block inherits. + +The fix, now implemented, is to make the aligner boundary #emph[canonical]: order and name blocks by +a hash of their consensus sequence, so that everything the aligner sees is a pure function of the +block content, independent of input order, tree order and identifier assignment. A second, initially +underestimated defect had to be fixed alongside it - neighbor joining resolves its Q-matrix ties by +matrix row order, which for three taxa is always tied, so the inferred topology itself followed the +argument order. All six permutations now agree, with or without a guide tree. + +Two claims in the first draft of this note did not survive measurement, and are corrected in place: +the mmseqs backend did #emph[not] share the defect, and neighbor-joining ties are not rare. The +first of those led to a simpler alternative --- enabling minimap2's dual mappings --- which was +tested, works, and is benchmarked against the landed fix in the final section. That benchmark in +turn corrected a third claim of this note's own: the apparent quality advantage of dual mappings was +a confounded comparison, and disappears once only one variable is changed. + +== Evidence + +=== The guide tree is obeyed + +Running orders `ERS990151 AP025961 NZ_CP035927` and `NZ_CP035927 ERS990151 AP025961` with the same +`--guide-tree` yields identical logged topology and an identical merge schedule: + +``` +Guide tree (newick): ((ERS990151,AP025961),NZ_CP035927); +=== Graph merging start: clades sizes 1 + 1 +=== Graph merging completed: clades sizes 1 + 1 -> 2 +=== Graph merging start: clades sizes 2 + 1 +=== Graph merging completed: clades sizes 2 + 1 -> 3 +``` + +`build_tree_from_newick` (`tree/newick.rs:70`) attaches graphs to leaves by name through a +`BTreeMap`, so topology and left/right assignment are fixed entirely by the tree file. Merge +scheduling is not the culprit. + +=== Order dependence survives when merge order cannot vary + +With two genomes and the tree `(ERS990151,AP025961)` there is exactly one possible merge. The +result still depends on argument order: + +#table( + columns: (auto, auto, auto, auto), + stroke: none, + table.header( + [*argument order*], [*blocks*], [*consensus shared with other order*], [*total consensus bp*], + ), + table.hline(), + [`ERS990151 AP025961`], [62], [37 / 62], [3 358 923], + [`AP025961 ERS990151`], [62], [37 / 62], [3 358 784], +) + +Same block count, same depth multiset, but 25 of 62 blocks carry a different consensus sequence and +the total consensus length differs by 139 bp. This rules out merge scheduling conclusively. + +=== The alignment direction flips + +A probe calling `align_with_minimap2_lib` directly on the two genomes, varying only the `BlockId` +assignment, gives: + +#table( + columns: (auto, auto, auto, auto), + stroke: none, + table.header( + [*ids*], [*hits*], [*cross-genome direction (qry #sym.arrow ref)*], [*matched bp*], + ), + table.hline(), + [0, 1], [217], [`ERS990151 -> AP025961` #sym.times 147], [263 139], + [1, 0], [216], [`AP025961 -> ERS990151` #sym.times 146], [266 047], + [1, 2], [217], [`ERS990151 -> AP025961` #sym.times 147], [263 139], + [2, 10], [216], [`AP025961 -> ERS990151` #sym.times 146], [266 047], +) + +Three things follow. The cross-genome hits flip direction wholesale. The hit sets genuinely differ +(147 vs 146 hits, #sym.tilde 2 900 matched bases), so this is not a relabelling. And only the +#emph[relative] order matters --- `0,1` behaves identically to `1,2`, while `2,10` behaves like +`1,0`, proving the comparison is lexicographic on the decimal string rather than numeric, since +`"10" < "2"`. + +== Root cause + +A four-step chain from argument position to alignment direction. + ++ *Identifiers are the input position.* `Pangraph::singleton` (`pangraph/pangraph.rs:29`) sets + `NodeId(fasta.index)`, `BlockId(fasta.index)` and `PathId(fasta.index)`, where `fasta.index` is + the record's position in the concatenated input. + ++ *Identifiers become aligner sequence names.* `align_with_minimap2_lib` + (`align/minimap2_lib/align_with_minimap2_lib.rs:19`) stringifies the `BlockId` as the name passed + to the index. + ++ *`minimap2` picks the direction by string comparison.* The call sets `X: true` + (same file, line 55), which raises `MM_F_NO_DUAL`. In the vendored source, `skip_seed` + (`packages/minimap2-sys/minimap2/map.c:89`) drops every hit where + `strcmp(qname, tname) > 0`; the header documents the flag as + #emph["skip pairs where query name is lexicographically larger than target name"] + (`minimap.h:11`). Each pair is aligned in one direction only. + ++ *Query and reference are not interchangeable.* `assign_anchor_block` + (`pangraph/reweave.rs:144`) resolves depth and ambiguity ties with `if ref_n <= qry_n`, so the + reference wins ties and its consensus becomes the anchor that the merged block inherits + (`MergePromise::solve_promise` returns `self.anchor_block`). + +This also explains the reported #emph["the file listed last dominates"] pattern: the last-listed +file receives the highest index, hence a lexicographically large name, hence it is usually the +reference, hence usually the anchor --- so its sequence dominates the merged consensus. + +=== A second, independent defect + +`align_with_minimap2_lib_impl` collects its results with `par_bridge()` (line 65). Rayon documents +that #emph["the resulting iterator is not guaranteed to keep the order of the original iterator"] +(`rayon-1.12.0/src/iter/par_bridge.rs:26`). `filter_matches` (`pangraph/graph_merging.rs:199`) then +applies a #emph[stable] sort by energy and greedily accepts non-overlapping matches, so equal-energy +ties are resolved by an order Rayon explicitly does not guarantee. The reported control run was +byte-identical, so this is not what was observed --- but it is latent nondeterminism in the same +code path and should be closed alongside. + +== Why this needs fixing + +- *It changes the biology.* Block boundaries are where recombination breakpoints are called + downstream. A 12% swing in block count driven by argument order propagates directly into + biological conclusions. + +- *It defeats the guide tree's purpose.* Users pass `--guide-tree` precisely to control merging and + obtain comparable results. The current behaviour silently ignores that intent for the + query/reference decision. + +- *It is non-monotonic and unpredictable.* Because the comparison is lexicographic on decimal + strings, behaviour changes discontinuously as block counts cross powers of ten. There is no + mental model a user could form to anticipate it. + +- *It is a latent-defect generator.* Opaque identifiers are steering an algorithmic decision. Any + future change to identifier assignment silently changes scientific output. + +- *`-X` discards real signal.* The two directions find measurably different homology (147 vs 146 + hits). Currently one is chosen arbitrarily. This is a quality concern in its own right, pursued + in the dual-mapping alternative section. + +== The fix + +Attack the point where the identifier leaks into the algorithm: the aligner boundary. Order and +name blocks by a hash of their consensus, so that the entire aligner input is a pure function of the +block content multiset. + +=== Canonical block naming + +```rust +use std::collections::{BTreeMap, HashMap}; +use std::hash::{Hash, Hasher}; +use twox_hash::XxHash64; + +/// Content-derived naming and ordering for the blocks handed to an aligner. +/// +/// Names are fixed-width so that byte-wise (`strcmp`) comparison agrees with the +/// numeric `(consensus_hash, block_id)` order that `canonical_order` uses. +pub struct BlockNames { + names: BTreeMap, + ids: HashMap, + keys: BTreeMap, + order: Vec, +} + +/// Order-independent key for a block: hash of its consensus, tie-broken by id. +fn consensus_key(id: BlockId, block: &PangraphBlock) -> (u64, BlockId) { + let mut hasher = XxHash64::with_seed(0); + block.consensus().hash(&mut hasher); + (hasher.finish(), id) +} + +impl BlockNames { + pub fn from_blocks(blocks: &BTreeMap) -> Self { + let keys: BTreeMap = + blocks.iter().map(|(&id, b)| (id, consensus_key(id, b))).collect(); + + let mut order: Vec = keys.keys().copied().collect(); + order.sort_by_key(|id| keys[id]); + + // 16 hex digits for the hash, 20 decimal digits for the id (usize::MAX has 20). + // Both zero-padded, so lexicographic order == numeric order at every position. + let names: BTreeMap = keys + .iter() + .map(|(&id, &(h, _))| (id, format!("{h:016x}_{:020}", id.0))) + .collect(); + + let ids = names.iter().map(|(&id, n)| (n.clone(), id)).collect(); + + Self { names, ids, keys, order } + } + + /// Blocks in canonical (content-derived) order. + pub fn canonical_order(&self) -> impl Iterator + '_ { + self.order.iter().copied() + } + + pub fn name(&self, id: BlockId) -> &str { + &self.names[&id] + } + + /// Recover the `BlockId` an aligner hit refers to. Fallible: an unknown name means + /// the aligner returned something we did not feed it. + pub fn id_of(&self, name: &str) -> Result { + self + .ids + .get(name) + .copied() + .ok_or_else(|| make_internal_report!("Aligner returned unknown sequence name '{name}'")) + } + + /// Order-independent sort key, for canonical tie-breaking downstream. + pub fn sort_key(&self, id: BlockId) -> (u64, BlockId) { + self.keys[&id] + } +} +``` + +=== Plugging it into the aligner + +```rust +pub fn align_with_minimap2_lib( + blocks: &BTreeMap, + names: &BlockNames, + params: &AlignmentArgs, +) -> Result, Report> { + // Canonical order, canonical names: nothing here depends on BlockId assignment. + let (seq_names, seqs): (Vec<&str>, Vec<&str>) = names + .canonical_order() + .map(|id| (names.name(id), blocks[&id].consensus().as_str())) + .unzip(); + + align_with_minimap2_lib_impl(&seqs, &seq_names, names, params) +} +``` + +=== Deterministic parallel output ordering + +Replace the unordered bridge with an indexed parallel iterator. `par_iter` over a slice is an +`IndexedParallelIterator`, for which `collect` into a `Vec` preserves input order by construction: + +```rust + let results: Vec = seqs + .par_iter() + .zip(names_v.par_iter()) + .map_init( + || Minimap2Mapper::new(&idx).unwrap(), + |mapper, (seq, name)| { + mapper + .run_map(seq, name) + .wrap_err_with(|| format!("When aligning sequence '{name}'")) + }, + ) + .collect::, Report>>()?; +``` + +This is not merely a determinism fix: indexed splitting also gives Rayon better work division than +`par_bridge`, which has to serialise pulls from the underlying sequential iterator behind a mutex. + +`self_merge`'s other parallel stage, `mergers.into_par_iter()` (`pangraph/graph_merging.rs:145`), +is already indexed (`Vec::into_par_iter`) and needs no change. + +=== Canonical downstream tie-breaks + +Two places consume the query/reference distinction and must stop privileging the reference. + +Energy sorting in `filter_matches`, currently a stable sort whose ties fall through to vector order: + +```rust + let alns = alns + .iter() + .map(|aln| (aln, alignment_energy2(aln, args))) + .filter(|(_, energy)| *energy < 0.0) + .sorted_by(|(a, ea), (b, eb)| { + OrderedFloat(*ea) + .cmp(&OrderedFloat(*eb)) + .then_with(|| names.sort_key(a.qry.name).cmp(&names.sort_key(b.qry.name))) + .then_with(|| a.qry.interval.start.cmp(&b.qry.interval.start)) + .then_with(|| names.sort_key(a.reff.name).cmp(&names.sort_key(b.reff.name))) + .then_with(|| a.reff.interval.start.cmp(&b.reff.interval.start)) + }) + .map(|(aln, _)| aln) + .collect_vec(); +``` + +Anchor selection in `assign_anchor_block`: + +```rust +fn assign_anchor_block(mergers: &mut [Alignment], graph: &Pangraph, names: &BlockNames) { + for m in mergers.iter_mut() { + let ref_block = &graph.blocks[&m.reff.name]; + let qry_block = &graph.blocks[&m.qry.name]; + + let n_of = |b: &PangraphBlock, iv: &Interval| { + b.consensus()[iv.to_range()].iter().filter(|c| c.0 == b'N').count() + }; + + // Deeper block wins; then fewer ambiguous bases; then a content-derived key. + // The final arm replaces `ref_n <= qry_n`, which privileged the reference and + // therefore leaked input order into the choice of surviving consensus. + let anchor = ref_block + .depth() + .cmp(&qry_block.depth()) + .then_with(|| n_of(qry_block, &m.qry.interval).cmp(&n_of(ref_block, &m.reff.interval))) + .then_with(|| names.sort_key(m.qry.name).cmp(&names.sort_key(m.reff.name))); + + m.anchor_block = Some(match anchor { + Ordering::Less => AnchorBlock::Qry, + _ => AnchorBlock::Ref, + }); + } +} +``` + +=== Why this is stronger than ordering identifiers by the guide tree + +An obvious alternative is to assign singleton identifiers in tree-leaf order rather than FASTA +order when `--guide-tree` is given. It does achieve argument-order independence, but it is +strictly weaker: + +- It only covers the `--guide-tree` path, leaving the default invocation broken. The + neighbor-joining path has its own, independent order dependence, addressed separately below. +- It trades one arbitrary convention for another. `((A,B),C)` and `((B,A),C)` are the same topology + but would still give different graphs. +- It has an implementation trap: `PathId.0` is used as a direct index into the FASTA vector + (`reconstruct/reconstruct_run.rs:62` together with `build/build_run.rs:44`), so renumbering + without permuting the records silently verifies against the wrong sequences. + +By contrast, canonicalising the aligner boundary removes the dependence on the first two at once. +Because `graph_join` is symmetric (`map_merge` over `BTreeMap`s) and `merge_graphs` uses left/right +only for debug logging, sibling order in the tree also stops mattering once the aligner input is +canonical --- confirmed by measurement. The neighbor-joining path needs its own fix, which is the +one place where this approach is not sufficient on its own. + +== Problems this could introduce + +=== Hash collisions between identical blocks + +This is the sharpest hazard, and naming purely by content hash would be a regression. + +Two distinct blocks can legitimately carry an identical consensus --- repeated elements, or +duplicated regions in the same genome that have not yet merged. If both were given the same name, +`skip_seed` (`map.c:81`) would take the `MM_F_NO_DIAG` branch: + +```c +cmp = strcmp(qname, s->name); +if ((flag&MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen) { + if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagonal anchors + ... +} +``` + +For two identical sequences the true alignment #emph[is] the diagonal, so every anchor supporting it +would be discarded and the pair would never merge --- precisely the pair we most want merged. +`find_matches` would compound this: it filters `m.qry.name != m.reff.name` to drop self-alignments, +which with colliding names would also drop the genuine cross-block hit. + +The `{hash}_{id}` suffix resolves this. Distinct blocks always have distinct names, so `cmp == 0` +occurs only for a true self-comparison, and `NO_DIAG` retains exactly its intended meaning. + +The residual is bounded and benign. When two blocks share a consensus, the direction falls back to +`BlockId`, which is still input-order dependent. But the two sequences are identical, so the +alignment is symmetric and the anchor choice affects only which #emph[identifier] survives, not the +consensus. Downstream, identifiers no longer influence alignment, so the effect is confined to +labels in the output JSON. This should be stated in the docs rather than papered over: invariance is +guaranteed when consensus sequences are distinct. + +A true 64-bit `XxHash64` collision between #emph[different] sequences is a separate matter. At +#sym.tilde $10^4$ blocks the birthday probability is #sym.tilde $10^(-11)$, and the consequence is +merely a different-but-still-deterministic ordering, not incorrect output, since the identifier +recovered from the name is exact. No mitigation needed. + +=== Consequences of the identifier lookup table + +Replacing `BlockId::from_str(&paf.q.name)` with a table lookup has a wider blast radius than it +first appears: + +- *Signature churn.* `Alignment::from_minimap_paf_obj` + (`align/minimap2_lib/align_with_minimap2_lib.rs:89`) and `Alignment::from_paf_str` + (`align/mmseqs/paf.rs:40`) must both take `&BlockNames`. `find_matches` and `filter_matches` + (`pangraph/graph_merging.rs:176`, `:187`) gain the parameter, and `self_merge` constructs the + table once per iteration from `graph.blocks`. + +- *The mmseqs path breaks silently otherwise.* `PafTsvRecord` deserialises `query: BlockId` and + `target: BlockId` directly through serde (`align/mmseqs/paf.rs:15`, `:20`), which only works + while names are bare decimal integers. With canonical names this must become a `String` field + parsed through `id_of`. Missing this would turn a working backend into a parse error --- caught by + compilation only if the field type is changed deliberately. + +- *Lookups become fallible.* An unknown name currently cannot happen; with a table it can, so the + error path needs a real internal error rather than an `unwrap`. This is a strict improvement in + diagnosability. + +- *Cost is negligible.* The table is `O(n_blocks)` entries (hundreds to low thousands), rebuilt once + per `self_merge` iteration. Hashing all consensuses costs one pass over the block sequence set + --- a few Mbp at `XxHash64` throughput, i.e. milliseconds against multi-second alignment stages. + +=== The mmseqs backend did *not* carry the defect + +An earlier draft of this note claimed that `align_with_mmseqs` carried the identical defect by a +different route, since it also writes its FASTA in `BTreeMap` order under `id.to_string()` names. +Measurement contradicts that. On the same three *Campylobacter* genomes that swing minimap2 from +508 to 462 blocks, the mmseqs backend was already order-invariant before any change: + +#table( + columns: (auto, auto, auto), + stroke: none, + table.header( + [*backend / argument order*], [*blocks (before fix)*], [*consensus-set digest*], + ), + table.hline(), + [mmseqs, `ERS AP NZ`], [297], [`865e8320ea1d72b1`], + [mmseqs, `NZ ERS AP`], [297], [`865e8320ea1d72b1`], +) + +The reason is precisely the `-X` flag. mmseqs does *not* pass it, so it computes both directions of +every pair and the resulting alignment set is symmetric - independent of which sequence is nominally +the query. `filter_matches` then ranks by energy, and the ranking is itself order-independent, so no +tie-break is ever reached. `MM_F_NO_DUAL` is not merely one mechanism among several: it is *the* +mechanism, and only the minimap2 path used it. + +The canonicalization is still applied to mmseqs, for two reasons: it removes the backend's residual +reliance on vector-order tie-breaking in `filter_matches` (latent rather than observed), and it keeps +the two backends behaving identically rather than accidentally-equivalently. But it should be +recorded as hardening, not as a bug fix. + +It is not free, though. Because the anchor tie-break changed, mmseqs output *moves*: 297 blocks +before, 317 after, order-invariant in both cases. That is a different-but-equally-valid anchor +choice on a backend that was not broken, and anyone comparing mmseqs graphs across this change +should expect the shift. + +This diagnosis has a direct consequence: if dual mapping is what protected mmseqs, then enabling it +for minimap2 should protect minimap2 too, without any canonical naming. That prediction was tested, +and it holds --- see the dual-mapping alternative section. + +=== Downstream fixtures and `pypangraph` + +`pypangraph` needs no code change. Its schema (`pypangraph/pangraph_schema.py`) is autogenerated +from the Rust types and covers only the serialised graph --- `Pangraph`, `PangraphPath`, +`PangraphBlock`, `PangraphNode`, `Edit` --- with no reference to `Alignment`, `Hit` or the PAF +types. `BlockId` remains a `uint`, so the loader, `IndexedCollection` and the integer/string +identifier duality are all untouched. + +What does need attention is #emph[fixture regeneration]. `packages/pypangraph/tests/data/plasmids.json` +is committed data, so it does not move when the build changes --- but the moment it is regenerated +with a fixed binary, block boundaries, block counts and identifiers all shift, and the hard-coded +expectations in `tests/test_graph.py` break: + +#table( + columns: (auto, auto), + stroke: none, + table.header( + [*location*], [*pinned expectation*], + ), + table.hline(), + [`test_graph.py:101`], [literal block id `"14710008249239879492"`], + [`test_graph.py:61-63`], [137 blocks, 27 core, 10 duplicated], + [`test_graph.py:88-90`], [137 #sym.times 15 matrix, sum 1042], +) + +The same caveat applies to `data/test_graph.json` on the Rust side if it is ever rebuilt rather than +merely consumed. Neither fixture is regenerated by this change, so nothing breaks on landing; the +point is that regeneration must be a deliberate, separate step, with the assertions above updated in +the same commit rather than rediscovered as a mysterious test failure later. + +Worth recording in `pypangraph`'s changelog regardless: even after the structural guarantee lands, +`BlockId` values still derive from `fasta.index`, so downstream code must not assume identifiers are +stable across runs with different input order. + +=== Things that are safe + +- *Correctness of merging is unaffected.* Which block is anchor changes which valid consensus is + retained, not whether the result is valid. Sequence reconstruction and the existing sanity checks + are indifferent to the choice. + +- *Integration tests do not exercise this path.* `data/test_graph.json` is consumed by the export + tests (`itest_export_*.rs`), which read a pre-built graph. They will not shift. + +- *No serialised format changes.* `Hit.name` remains a `BlockId`; canonical names exist only for the + duration of an alignment call and never reach disk. + +- *No performance regression.* `-X` still halves the pair count; the index size and mapping work are + unchanged; the parallel stage gets marginally better splitting. + +=== Neighbor-joining ties are not rare, and had to be fixed too + +An earlier draft listed the neighbor-joining tie-break as an out-of-scope residual, on the +assumption that exact ties in the Q matrix are rare with real-valued mash distances. That reasoning +was wrong, and for the commonest small case it is wrong by construction. + +For three taxa, write $S_k$ for the sum of row $k$ of the distance matrix. With $n = 3$ the code +computes $Q_(i j) = D_(i j) - S_i - S_j$, so + +$ Q_12 = D_12 - (D_12 + D_13) - (D_12 + D_23) = -(D_12 + D_13 + D_23) $ + +and the same value comes out for $Q_13$ and $Q_23$. #emph[Every] off-diagonal entry is identical, +regardless of the data. `Q.argmin()` therefore always returns `(0, 1)`, and the tree is always +`((first, second), third)` in argument order. This is exactly what the reproducer shows: the same +three genomes give `((ERS990151,AP025961),NZ_CP035927)` in one order and +`((NZ_CP035927,ERS990151),AP025961)` in another - genuinely different topologies, not sibling swaps. + +Since a different topology legitimately yields a different graph, canonicalizing the aligner +boundary cannot help here; the tree itself has to stop depending on argument order. The fix is to +order the leaves handed to neighbor joining by a content-derived key, so that the row order the +tie-break falls back on is a property of the sequences rather than of the command line: + +```rust +fn graph_content_key(graph: &Pangraph) -> Vec { + graph.blocks.values().map(consensus_hash).sorted().collect() +} + +pub fn build_tree_using_neighbor_joining(graphs: Vec) -> Result<...> { + let mut graphs = graphs; + graphs.sort_by_cached_key(graph_content_key); + // ... +} +``` + +Without this, the default invocation - no `--guide-tree` - remains order-dependent, so a fix that +stopped at the aligner boundary would have left the common case broken. + +=== Residual order dependence after the fix + +Honesty about what is #emph[not] fixed: + ++ Identical-consensus blocks, as discussed, retain identifier-level (not sequence-level) order + dependence. ++ `BlockId` values themselves still derive from `fasta.index`, so the identifiers appearing in + output JSON differ between argument orders even when the graph structure is identical. If + byte-identical output is wanted, singleton identifiers would also need a content-derived + assignment --- worth considering, but orthogonal to the correctness issue. + +With those caveats, the guarantee becomes: #emph[the graph structure is a pure function of the input +sequence set, the tree topology, and the alignment parameters] --- independent of argument order and +of sibling order within the tree. + +== Validation results + +Implemented on branch `debug/order-dependence` across five commits: `BlockNames` and its unit tests, +the aligner boundary, the downstream tie-breaks, the neighbor-joining leaf ordering, and the +invariance regression tests. + +=== The reproducer collapses + +All six permutations of the three *Campylobacter* genomes, `--circular --guide-tree`: + +#table( + columns: (auto, auto, auto), + stroke: none, + table.header( + [*argument order*], [*blocks before*], [*blocks after*], + ), + table.hline(), + [`ERS AP NZ`], [508], [512], + [`ERS NZ AP`], [462], [512], + [`AP ERS NZ`], [497], [512], + [`AP NZ ERS`], [452], [512], + [`NZ ERS AP`], [462], [512], + [`NZ AP ERS`], [452], [512], +) + +The consensus multisets are byte-identical across all six (SHA-256 `5495a5963cd3b2b8`), as are the +depth multisets. A sibling-swapped guide tree gives the same digest, and the neighbor-joining path +(no `--guide-tree`) likewise collapses to a single answer across all six orders - 522 blocks, which +differs from 512 only because it infers a different topology, as it should. + +The JSON is *not* byte-identical between orders, exactly as predicted: `BlockId` still derives from +`fasta.index`, so identifiers move even when structure does not. + +=== The regression tests discriminate + +`packages/pangraph/tests/itest_build_order_invariance.rs` asserts invariance over six argument +permutations under three topologies, plus sibling-swap invariance. Each assertion was checked +against deliberately reverted code: + +#table( + columns: (auto, auto), + stroke: none, + table.header( + [*reverted component*], [*failing case*], + ), + table.hline(), + [canonical names + anchor tie-break], [`balanced_guide_tree`, order `[3,2,1,0]`], + [neighbor-joining leaf ordering], [`no_guide_tree`, order `[3,2,1,0]`], +) + +One lesson is worth recording, because the first version of this test was worthless. A +#emph[three]-genome fixture cannot detect the anchor defect: the final merge joins a depth-2 clade +with a depth-1 leaf, so depth decides and the tie-break never runs. The test passed with the fix +fully reverted. It needs four genomes and a balanced topology, so that the top merge is depth-2 +against depth-2 - a genuine tie. Any future fixture for this class of bug must be checked against +reverted code before it is trusted. + +=== Fixtures are untouched + +`data/test_graph.json` and `packages/pypangraph/tests/data/plasmids.json` are byte-identical after +the change, so both suites verify the fix rather than absorb it. All 327 unit tests and all +integration tests pass; `clippy --all-targets -Dwarnings` is clean. + +=== Still unverified + +The two-identical-consensus-blocks merge case is covered only at the naming level +(`identical_consensus_blocks_get_distinct_names`), not end to end through a real alignment. The +collision hazard argument in the corresponding section is therefore reasoned, not measured. + +== The dual-mapping alternative + +The finding that mmseqs was never affected raises an obvious question: if computing both directions +is what made mmseqs immune, could minimap2 simply be told to do the same --- fixing the defect +without any of the naming machinery? The answer is yes, with caveats that decide whether it is the +right trade. + +=== Disabling `-X` outright does not work + +`X` is a bundle, not a toggle. In the in-tree wrapper (`packages/minimap2/src/options_args.rs:324`) +a single boolean sets four flags: + +```rust + if args.X { + map_opt.flag |= (MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN) as i64; + } +``` + +and the `asm5`/`asm10`/`asm20` presets pangraph uses set none of them (unlike the `ava-*` presets, +`packages/minimap2-sys/minimap2/options.c:96`). So `X: false` also discards all-chains and +no-diagonal reporting. Measured on the three *Campylobacter* genomes, the result is order-invariant +and useless: + +#table( + columns: (auto, auto, auto), + stroke: none, + table.header( + [*variant*], [*blocks*], [*invariant across six orders*], + ), + table.hline(), + [`X: false`], [3], [yes], + [`-X` with only `MM_F_NO_DUAL` cleared], [459], [yes], +) + +Three blocks for three genomes means nothing merged at all. The likely mechanism --- inferred from +minimap2's secondary-alignment filtering rather than instrumented --- is that without `NO_DIAG` each +sequence's perfect self-hit becomes the primary chain, and cross-genome hits are then discarded as +low-ratio secondaries. + +=== Enabling dual mappings does work + +The precise variant is minimap2's `-X --dual=yes`: keep all-chains, no-diagonal and no-long-join, +clear only `MM_F_NO_DUAL`. Tested against the #emph[entire pre-fix codebase] --- `BlockId` sequence +names, `BlockId` ordering, reference-wins-ties anchor selection, input-order neighbor-joining leaves +--- all six argument permutations produced 459 blocks with an identical consensus-set digest +(`e33506c21d945b73`). No canonical naming, no lookup table, no tie-break changes. + +This is a genuinely simpler fix. It needs one wrapper change: split `X: bool` into `X` plus a +separate `dual` option, since nothing currently exposes `--dual`. + +=== A confounded comparison, corrected + +An earlier version of this section reported that dual mappings gave 459 blocks against 512 for the +landed fix, and reasoned from that gap that dual mappings produce a better-consolidated graph. That +comparison was #emph[confounded]: it set pre-fix code #sym.plus dual mappings against canonical +naming #sym.plus `-X`, changing two things at once. The 459 was a property of the pre-fix anchor +tie-break, not of dual mappings. + +Holding canonical naming fixed and toggling only dual mappings, the same three genomes give 512 +blocks against 514 --- a difference of two blocks, not fifty-three. The quality argument built on +the larger gap does not survive. + +=== Benchmark across the bundled datasets + +Both variants were run over every dataset in `data/`, plus the three-genome reproducer, with +canonical naming held fixed so that dual mapping is the only variable. Graph metrics are computed +from the JSON directly; the definitions were cross-checked against +`pypangraph.to_blockstats_df()` on `data/test_graph.json` and agree exactly (14 blocks, 6 core, +33 611 core bp, 50 791 pangenome bp). A block counts as #emph[core] when every path crosses it +exactly once. + +#figure( + image("assets/n00/dual_quality.png", width: 92%), + caption: [Graph metrics are unchanged on seven of ten datasets --- `example`, `flu-h1`, `flu-h3`, + `ges-1`, `mpox`, `russian_doll_plasmids` and `sc2` are bit-identical. Where they do move, every + change is under 1%, and not consistently signed.], +) + +#figure( + image("assets/n00/dual_runtime.png", width: 92%), + caption: [Runtime overhead of dual mappings. Datasets completing in under a second are omitted as + timing noise. Measured with nothing else running on the machine.], +) + +Absolute values for the datasets where anything changed: + +#table( + columns: (auto, auto, auto, auto, auto), + stroke: none, + table.header( + [*dataset*], [*variant*], [*blocks*], [*core bp*], [*pangenome bp*], + ), + table.hline(), + [klebs (9 genomes)], [`-X`], [1 381], [4 458 249], [7 644 985], + [], [dual], [1 372], [4 457 866], [7 644 453], + [ecoli (10 genomes)], [`-X`], [2 914], [3 782 006], [7 830 290], + [], [dual], [2 908], [3 782 120], [7 827 250], + [campylobacter (3)], [`-X`], [512], [102 659], [3 445 626], + [], [dual], [514], [102 867], [3 445 993], +) + +Two things follow. Dual mappings are close to free in graph terms --- seven of ten datasets are +bit-identical, and the three that move do so by well under 1%. And the direction of the change is +inconsistent: dual mappings #emph[reduce] the block count on klebs (#sym.minus 0.65%) and ecoli +(#sym.minus 0.21%) but #emph[increase] it on Campylobacter (#sym.plus 0.39%). There is no evidence +here that either variant consolidates the graph better; they simply make marginally different +choices at a handful of boundaries. + +The cost, by contrast, is real and consistent: every dataset gets slower, by 5% (`ges-1`) to 26% +(Campylobacter), with the two large bacterial pangenomes at #sym.plus 11% and #sym.plus 19%. That is +well short of the 2#sym.times one might assume --- alignment is not the whole pipeline --- but it is +a uniform tax for no measured graph benefit. + +The remaining argument for dual mappings is therefore the one from evidence rather than from +outcome: the direction probe showed the two directions genuinely find different homology (147 versus +146 hits, 263 139 versus 266 047 matched bases), and under `-X` one of them is discarded +arbitrarily. Using both is more principled. It just does not, on these datasets, produce a +measurably different graph. + +=== What it does not buy + +The decisive distinction is elsewhere. Dual mappings do not remove the order-sensitive tie-breaks; +they merely stop them being reached. The energy sort still falls through to vector order on exact +ties, and reference-wins-ties would still privilege one side. Mirror alignments happen to score +slightly differently, so no tie arises in practice --- which is exactly why mmseqs looked immune. +#emph[Dual mappings make invariance a property of the data; canonical naming makes it a property of +the code.] If two mirror alignments ever scored equal, the order dependence would return. + +=== Recommendation + +Keep the canonical naming, and do not adopt dual mappings on the strength of this evidence. + +The case for dual mappings rested on two claims, and the benchmark removes both. It was to be the +simpler fix --- but it is not sufficient on its own, since the neighbor-joining ordering still has +to be fixed separately, and it only makes invariance contingent on mirror alignments never tying. +And it was to produce a better graph --- but across ten datasets it produces the #emph[same] graph +seven times, sub-1% differences three times, in no consistent direction, at a uniform 5--26% runtime +cost. + +What remains is a principled objection to `-X`: discarding one direction of a pairwise alignment +throws away real signal. That objection stands, and if the alignment stage is ever revisited for +quality reasons, dual mappings are the right thing to reach for --- ideally with a quality metric +better than block counts to judge by. On present evidence it is a cost with no measured benefit, +and the canonical naming already delivers the guarantee it was meant to provide. diff --git a/packages/pangraph/src/align/block_names.rs b/packages/pangraph/src/align/block_names.rs new file mode 100644 index 00000000..d490c888 --- /dev/null +++ b/packages/pangraph/src/align/block_names.rs @@ -0,0 +1,203 @@ +use crate::make_internal_report; +use crate::pangraph::pangraph_block::{BlockId, PangraphBlock}; +use eyre::Report; +use std::collections::{BTreeMap, HashMap}; +use std::hash::{Hash, Hasher}; +use twox_hash::XxHash64; + +/// Number of hex digits used to render the consensus hash in a canonical name. +const HASH_WIDTH: usize = 16; + +/// Number of decimal digits used to render the block id in a canonical name. +/// `usize::MAX` is 20 digits, so this never truncates. +const ID_WIDTH: usize = 20; + +/// Content-derived ordering key for a block: hash of its consensus sequence, tie-broken by block id. +/// +/// The tie-break only engages for blocks whose consensus is byte-identical, in which case the +/// alignment between them is symmetric and the choice cannot bias the result. +pub type BlockKey = (u64, BlockId); + +/// Content-derived naming and ordering for the blocks handed to an alignment backend. +/// +/// `BlockId`s are assigned from the position of a record in the input (see `Pangraph::singleton`), +/// so using them to name sequences leaks the input order into the aligner. `minimap2` run with +/// `-X` keeps only one direction of each pair, chosen by `strcmp` of the sequence names, which +/// makes the query/reference roles - and therefore the merged consensus - depend on the order the +/// FASTA files were listed on the command line. +/// +/// This type replaces those names with a key derived from the block consensus, so that both the +/// order blocks are fed to the aligner and the names they carry are a pure function of the block +/// content. Names are zero-padded to a fixed width, so byte-wise (`strcmp`) comparison agrees with +/// the numeric `(consensus_hash, block_id)` order that [`BlockNames::canonical_order`] uses. +pub struct BlockNames { + names: BTreeMap, + ids: HashMap, + keys: BTreeMap, + order: Vec, +} + +/// Hashes a block's consensus sequence. +/// +/// Content-derived, and therefore independent of how `BlockId`s were assigned. +pub fn consensus_hash(block: &PangraphBlock) -> u64 { + let mut hasher = XxHash64::with_seed(0); + block.consensus().hash(&mut hasher); + hasher.finish() +} + +/// Computes the content-derived ordering key of a block. +fn consensus_key(id: BlockId, block: &PangraphBlock) -> BlockKey { + (consensus_hash(block), id) +} + +/// Renders a canonical, fixed-width name for a block. +fn canonical_name((hash, id): BlockKey) -> String { + format!("{hash:0HASH_WIDTH$x}_{:0ID_WIDTH$}", id.0) +} + +impl BlockNames { + /// Builds the canonical naming and ordering for a set of blocks. + pub fn from_blocks(blocks: &BTreeMap) -> Self { + let keys: BTreeMap = blocks + .iter() + .map(|(&id, block)| (id, consensus_key(id, block))) + .collect(); + + let mut order: Vec = keys.keys().copied().collect(); + order.sort_unstable_by_key(|id| keys[id]); + + let names: BTreeMap = keys.iter().map(|(&id, &key)| (id, canonical_name(key))).collect(); + let ids: HashMap = names.iter().map(|(&id, name)| (name.clone(), id)).collect(); + + Self { + names, + ids, + keys, + order, + } + } + + /// Blocks in canonical (content-derived) order. + pub fn canonical_order(&self) -> impl Iterator + '_ { + self.order.iter().copied() + } + + /// Canonical name of a block, as handed to the alignment backend. + pub fn name(&self, id: BlockId) -> &str { + &self.names[&id] + } + + /// Recovers the block a canonical name refers to. + /// + /// Fails if the aligner returned a name we never fed it, which would otherwise surface as a + /// silently wrong block id. + pub fn id_of(&self, name: &str) -> Result { + self + .ids + .get(name) + .copied() + .ok_or_else(|| make_internal_report!("Aligner returned unknown sequence name '{name}'")) + } + + /// Content-derived sort key of a block, for order-independent tie-breaking. + pub fn sort_key(&self, id: BlockId) -> BlockKey { + self.keys[&id] + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::pangraph::pangraph_node::NodeId; + use crate::representation::seq::Seq; + use itertools::Itertools; + use pretty_assertions::assert_eq; + use rstest::rstest; + + /// Builds a block map from `(id, consensus)` pairs. + fn blocks_of(spec: &[(usize, &str)]) -> BTreeMap { + spec + .iter() + .map(|&(id, seq)| { + let bid = BlockId(id); + (bid, PangraphBlock::from_consensus(Seq::from_str(seq), bid, NodeId(id))) + }) + .collect() + } + + #[rstest] + fn canonical_names_are_fixed_width_and_unique() { + let blocks = blocks_of(&[(0, "ACGT"), (1, "TTTT"), (12345, "GGGG")]); + let names = BlockNames::from_blocks(&blocks); + + let rendered = blocks.keys().map(|&id| names.name(id).to_owned()).collect_vec(); + assert!(rendered.iter().all(|n| n.len() == HASH_WIDTH + 1 + ID_WIDTH)); + assert_eq!(rendered.iter().unique().count(), rendered.len()); + } + + #[rstest] + fn canonical_order_matches_lexicographic_name_order() { + let blocks = blocks_of(&[(0, "ACGT"), (1, "TTTT"), (2, "GGGG"), (3, "CCCC")]); + let names = BlockNames::from_blocks(&blocks); + + let by_order = names.canonical_order().map(|id| names.name(id)).collect_vec(); + let sorted = by_order.iter().copied().sorted().collect_vec(); + assert_eq!(by_order, sorted); + } + + /// The property the fix exists for: permuting block ids over the same consensus set must not + /// change the order or the relative naming the aligner sees. + #[rstest] + fn canonical_order_is_invariant_under_id_permutation() { + let seqs = ["ACGTACGT", "TTTTGGGG", "GGGGCCCC", "CACACACA"]; + + let forward = blocks_of(&seqs.iter().enumerate().map(|(i, s)| (i, *s)).collect_vec()); + let reversed = blocks_of( + &seqs + .iter() + .enumerate() + .map(|(i, s)| (seqs.len() - 1 - i, *s)) + .collect_vec(), + ); + + let n_fwd = BlockNames::from_blocks(&forward); + let n_rev = BlockNames::from_blocks(&reversed); + + // The sequences appear in the same order regardless of how ids were assigned. + let seq_of = |blocks: &BTreeMap, id: BlockId| blocks[&id].consensus().as_str().to_owned(); + let order_fwd = n_fwd.canonical_order().map(|id| seq_of(&forward, id)).collect_vec(); + let order_rev = n_rev.canonical_order().map(|id| seq_of(&reversed, id)).collect_vec(); + assert_eq!(order_fwd, order_rev); + } + + #[rstest] + fn names_round_trip_to_block_ids() { + let blocks = blocks_of(&[(0, "ACGT"), (7, "TTTT")]); + let names = BlockNames::from_blocks(&blocks); + + for &id in blocks.keys() { + assert_eq!(names.id_of(names.name(id)).unwrap(), id); + } + } + + #[rstest] + fn unknown_name_is_rejected() { + let blocks = blocks_of(&[(0, "ACGT")]); + let names = BlockNames::from_blocks(&blocks); + let err = names.id_of("not-a-name").unwrap_err(); + assert!(err.to_string().contains("unknown sequence name")); + } + + /// Blocks sharing a consensus must still receive distinct names: identical names would make + /// `minimap2` treat the pair as a self-comparison and discard the diagonal anchors that carry + /// their (perfect) alignment, so the two blocks could never merge. + #[rstest] + fn identical_consensus_blocks_get_distinct_names() { + let blocks = blocks_of(&[(0, "ACGTACGT"), (1, "ACGTACGT")]); + let names = BlockNames::from_blocks(&blocks); + + assert_ne!(names.name(BlockId(0)), names.name(BlockId(1))); + assert_eq!(names.sort_key(BlockId(0)).0, names.sort_key(BlockId(1)).0); + } +} diff --git a/packages/pangraph/src/align/minimap2_lib/align_with_minimap2_lib.rs b/packages/pangraph/src/align/minimap2_lib/align_with_minimap2_lib.rs index 348e4fb9..2db5c8bc 100644 --- a/packages/pangraph/src/align/minimap2_lib/align_with_minimap2_lib.rs +++ b/packages/pangraph/src/align/minimap2_lib/align_with_minimap2_lib.rs @@ -1,10 +1,11 @@ use crate::align::alignment::{Alignment, Hit}; use crate::align::alignment_args::AlignmentArgs; +use crate::align::block_names::BlockNames; use crate::pangraph::pangraph_block::{BlockId, PangraphBlock}; use crate::pangraph::strand::Strand; use crate::{make_error, make_internal_error}; use eyre::{Report, WrapErr}; -use itertools::{Itertools, izip}; +use itertools::Itertools; use minimap2::{Minimap2Args, Minimap2Index, Minimap2Mapper, Minimap2Preset, Minimap2Result}; use noodles::sam::record::Cigar; use num_traits::clamp_min; @@ -12,23 +13,33 @@ use rayon::prelude::*; use std::collections::BTreeMap; use std::str::FromStr; +/// Aligns the consensus sequences of `blocks` against each other, all-vs-all. +/// +/// Blocks are fed in canonical (content-derived) order under canonical names, so that neither the +/// set of alignments nor the query/reference role of each pair depends on how `BlockId`s were +/// assigned - and therefore on the order the input FASTA files were listed. See [`BlockNames`]. pub fn align_with_minimap2_lib( blocks: &BTreeMap, + names: &BlockNames, params: &AlignmentArgs, ) -> Result, Report> { - let (names, seqs): (Vec, Vec<&str>) = blocks - .iter() - .map(|(id, block)| (id.to_string(), block.consensus().as_str())) + let (seq_names, seqs): (Vec<&str>, Vec<&str>) = names + .canonical_order() + .map(|id| (names.name(id), blocks[&id].consensus().as_str())) .unzip(); - let alns: Vec = align_with_minimap2_lib_impl(&seqs, &names, params)?; + let alns: Vec = align_with_minimap2_lib_impl(&seqs, &seq_names, &|n| names.id_of(n), params)?; Ok(alns) } +/// Aligner mechanics, decoupled from how sequence names map back to blocks. +/// +/// `resolve` recovers the [`BlockId`] a name refers to; the caller decides the naming scheme. fn align_with_minimap2_lib_impl( seqs: &[impl AsRef], names: &[impl AsRef], + resolve: &(dyn Fn(&str) -> Result + Sync), params: &AlignmentArgs, ) -> Result, Report> { if names.len() != seqs.len() { @@ -61,11 +72,15 @@ fn align_with_minimap2_lib_impl( let idx = Minimap2Index::new(&seqs, &names, &args)?; - let results: Vec = izip!(&seqs, &names) - .par_bridge() + // `par_iter().zip()` over slices is an indexed parallel iterator, for which `collect` preserves + // input order by construction. `par_bridge()` does not guarantee order, which would leave the + // energy-sort tie-breaking in `filter_matches` at the mercy of thread scheduling. + let results: Vec = seqs + .par_iter() + .zip(names.par_iter()) .map_init( || Minimap2Mapper::new(&idx).unwrap(), - move |mapper, (seq, name)| { + |mapper, (seq, name)| { mapper .run_map(seq, name) .wrap_err_with(|| format!("When aligning sequence '{name}'")) @@ -75,7 +90,7 @@ fn align_with_minimap2_lib_impl( let alns = results .into_iter() - .map(Alignment::from_minimap_paf_obj) + .map(|res| Alignment::from_minimap_paf_obj(res, resolve)) .collect::>, Report>>()? .into_iter() .flatten() @@ -86,7 +101,10 @@ fn align_with_minimap2_lib_impl( #[allow(clippy::multiple_inherent_impl)] impl Alignment { - pub fn from_minimap_paf_obj(res: Minimap2Result) -> Result, Report> { + pub fn from_minimap_paf_obj( + res: Minimap2Result, + resolve: &(dyn Fn(&str) -> Result + Sync), + ) -> Result, Report> { let Minimap2Result { pafs, .. } = res; pafs .into_iter() @@ -94,12 +112,12 @@ impl Alignment { if let Some(cg) = &paf.cg { Ok(Alignment { qry: Hit::new( - BlockId::from_str(&paf.q.name)?, + resolve(&paf.q.name)?, paf.q.len, (paf.q.start as usize, paf.q.end as usize), ), reff: Hit::new( - BlockId::from_str(&paf.t.name)?, + resolve(&paf.t.name)?, paf.t.len, (paf.t.start as usize, paf.t.end as usize), ), @@ -183,7 +201,8 @@ mod tests { ..AlignmentArgs::default() }; - let actual = align_with_minimap2_lib_impl(&seqs, &names, ¶ms)?; + // Names here are the FASTA record ids, so resolve them as plain integers. + let actual = align_with_minimap2_lib_impl(&seqs, &names, &|n| BlockId::from_str(&n.to_owned()), ¶ms)?; let expected = vec![Alignment { qry: Hit::new(BlockId(0), 998, (0, 996)), diff --git a/packages/pangraph/src/align/mmseqs/align_with_mmseqs.rs b/packages/pangraph/src/align/mmseqs/align_with_mmseqs.rs index 132ae5ea..434801cb 100644 --- a/packages/pangraph/src/align/mmseqs/align_with_mmseqs.rs +++ b/packages/pangraph/src/align/mmseqs/align_with_mmseqs.rs @@ -1,5 +1,6 @@ use crate::align::alignment::Alignment; use crate::align::alignment_args::AlignmentArgs; +use crate::align::block_names::BlockNames; use crate::align::mmseqs::paf::PafTsvRecord; use crate::io::fasta::FastaWriter; use crate::io::file::open_file_or_stdin; @@ -15,8 +16,14 @@ use std::io::Read; use std::process::Command; use tempfile::Builder as TempDirBuilder; +/// Aligns the consensus sequences of `blocks` against each other using mmseqs. +/// +/// Like the minimap2 backend, blocks are written in canonical (content-derived) order under +/// canonical names, so the result does not depend on how `BlockId`s were assigned. See +/// [`BlockNames`]. pub fn align_with_mmseqs( blocks: &BTreeMap, + names: &BlockNames, params: &AlignmentArgs, ) -> Result, Report> { // TODO: This uses a global resource - filesystem. @@ -29,9 +36,9 @@ pub fn align_with_mmseqs( { let mut writer = FastaWriter::from_path(&input_path)?; - blocks - .iter() - .try_for_each(|(id, block)| writer.write(id.to_string(), &None, block.consensus()))?; + names + .canonical_order() + .try_for_each(|id| writer.write(names.name(id), &None, blocks[&id].consensus()))?; } let output_column_names = PafTsvRecord::fields_names().join(","); @@ -66,7 +73,7 @@ pub fn align_with_mmseqs( let mut paf_str = String::new(); open_file_or_stdin(&Some(output_path))?.read_to_string(&mut paf_str)?; - Alignment::from_paf_str(&paf_str) + Alignment::from_paf_str(&paf_str, &|n| names.id_of(n)) } // FIXME: This test is failing after commit a62b19b018b4b2f9602bc75335d4ab5ddbc7abf5 diff --git a/packages/pangraph/src/align/mmseqs/paf.rs b/packages/pangraph/src/align/mmseqs/paf.rs index 72af121a..c1c62162 100644 --- a/packages/pangraph/src/align/mmseqs/paf.rs +++ b/packages/pangraph/src/align/mmseqs/paf.rs @@ -12,12 +12,15 @@ use std::io::Cursor; #[allow(dead_code)] #[derive(Clone, Debug, Deserialize)] pub struct PafTsvRecord { - /* 01 */ query: BlockId, + // `query` and `target` are sequence *names*, resolved to `BlockId`s by the caller. They cannot + // be deserialized as `BlockId` directly, since canonical names are not bare integers. + /* 01 */ + query: String, /* 02 */ qlen: usize, /* 03 */ qstart: usize, /* 04 */ qend: usize, /* 05 */ empty: String, - /* 06 */ target: BlockId, + /* 06 */ target: String, /* 07 */ tlen: usize, /* 08 */ tstart: usize, /* 09 */ tend: usize, @@ -37,7 +40,11 @@ impl PafTsvRecord { #[allow(clippy::multiple_inherent_impl)] impl Alignment { - pub fn from_paf_str(paf_str: impl AsRef) -> Result, Report> { + /// Parses mmseqs PAF output, mapping sequence names back to blocks through `resolve`. + pub fn from_paf_str( + paf_str: impl AsRef, + resolve: &(dyn Fn(&str) -> Result + Sync), + ) -> Result, Report> { let mut rdr = CsvReaderBuilder::new() .delimiter(b'\t') .has_headers(false) @@ -53,8 +60,8 @@ impl Alignment { let (tstart, tend, _) = order_range(paf.tstart, paf.tend); Ok(Alignment { - qry: Hit::new(paf.query, paf.qlen, (qstart, qend)), - reff: Hit::new(paf.target, paf.tlen, (tstart, tend)), + qry: Hit::new(resolve(&paf.query)?, paf.qlen, (qstart, qend)), + reff: Hit::new(resolve(&paf.target)?, paf.tlen, (tstart, tend)), matches: paf.nident, length: paf.alnlen, quality: paf.bits, @@ -107,7 +114,10 @@ mod tests { divergence: Some(0.134), align: Some(693.0), }]; - assert_eq!(Alignment::from_paf_str(paf_content).unwrap(), aln); + assert_eq!( + Alignment::from_paf_str(paf_content, &|n| BlockId::from_str(&n.to_owned())).unwrap(), + aln + ); } #[rstest] @@ -127,6 +137,9 @@ mod tests { divergence: Some(0.134), align: Some(693.0), }]; - assert_eq!(Alignment::from_paf_str(paf_content).unwrap(), aln); + assert_eq!( + Alignment::from_paf_str(paf_content, &|n| BlockId::from_str(&n.to_owned())).unwrap(), + aln + ); } } diff --git a/packages/pangraph/src/align/mod.rs b/packages/pangraph/src/align/mod.rs index 8c7077e4..f5dda6bd 100644 --- a/packages/pangraph/src/align/mod.rs +++ b/packages/pangraph/src/align/mod.rs @@ -1,6 +1,7 @@ pub mod alignment; pub mod alignment_args; pub mod bam; +pub mod block_names; pub mod energy; pub mod map_variations; pub mod minimap2_lib; diff --git a/packages/pangraph/src/bin/mmseqs_example.rs b/packages/pangraph/src/bin/mmseqs_example.rs index f48ff4ff..13050f76 100644 --- a/packages/pangraph/src/bin/mmseqs_example.rs +++ b/packages/pangraph/src/bin/mmseqs_example.rs @@ -3,6 +3,7 @@ use ctor::ctor; use eyre::Report; use maplit::btreemap; use pangraph::align::alignment_args::AlignmentArgs; +use pangraph::align::block_names::BlockNames; use pangraph::align::mmseqs::align_with_mmseqs::align_with_mmseqs; use pangraph::io::fasta::FastaReader; use pangraph::pangraph::pangraph_block::{BlockId, PangraphBlock}; @@ -42,7 +43,8 @@ fn main() -> Result<(), Report> { .map(|block| (block.id(), block)) .collect(); - let result = align_with_mmseqs(&blocks, ¶ms)?; + let names = BlockNames::from_blocks(&blocks); + let result = align_with_mmseqs(&blocks, &names, ¶ms)?; println!("{:#?}", &result); Ok(()) diff --git a/packages/pangraph/src/pangraph/graph_merging.rs b/packages/pangraph/src/pangraph/graph_merging.rs index 5c2ad4a8..fad484cc 100644 --- a/packages/pangraph/src/pangraph/graph_merging.rs +++ b/packages/pangraph/src/pangraph/graph_merging.rs @@ -1,5 +1,6 @@ use crate::align::alignment::Alignment; use crate::align::alignment_args::AlignmentArgs; +use crate::align::block_names::BlockNames; use crate::align::energy::alignment_energy2; use crate::align::minimap2_lib::align_with_minimap2_lib::align_with_minimap2_lib; use crate::align::mmseqs::align_with_mmseqs::align_with_mmseqs; @@ -93,9 +94,14 @@ pub fn graph_join(left_graph: &Pangraph, right_graph: &Pangraph) -> Pangraph { } pub fn self_merge(graph: Pangraph, args: &PangraphBuildArgs) -> Result<(Pangraph, bool), Report> { + // Canonical, content-derived names and ordering for this round of alignment. Everything the + // aligner sees - and every tie-break that consumes its output - is keyed on these rather than + // on `BlockId`s, which carry the input order. + let names = BlockNames::from_blocks(&graph.blocks); + // use minimap2 or other aligners to find matches between the consensus // sequences of the blocks - let matches = find_matches(&graph.blocks, args)?; + let matches = find_matches(&graph.blocks, &names, args)?; debug!("Found matches: {}", matches.len()); trace!("{matches:#?}"); @@ -118,7 +124,7 @@ pub fn self_merge(graph: Pangraph, args: &PangraphBuildArgs) -> Result<(Pangraph // - calculate energy and keep only matches with E < 0 // - sort them by energy // - discard incompatible matches (the ones that have overlapping regions) - let mut matches = filter_matches(&matches, &args.aln_args); + let mut matches = filter_matches(&matches, &names, &args.aln_args); debug!("Matches after filtering: {}", matches.len()); trace!("{matches:#?}"); @@ -139,7 +145,7 @@ pub fn self_merge(graph: Pangraph, args: &PangraphBuildArgs) -> Result<(Pangraph // - adding the blocks that do not need merging to the preliminary graph // - return the set of blocks that should be merged let (mut graph, mergers) = - reweave(&mut matches, graph, args.aln_args.indel_len_threshold).wrap_err("During reweave")?; + reweave(&mut matches, graph, &names, args.aln_args.indel_len_threshold).wrap_err("During reweave")?; let mut merged_blocks: Vec = mergers .into_par_iter() @@ -175,16 +181,17 @@ pub fn self_merge(graph: Pangraph, args: &PangraphBuildArgs) -> Result<(Pangraph // Returns a list of alignment objects. pub fn find_matches( blocks: &BTreeMap, + names: &BlockNames, args: &PangraphBuildArgs, ) -> Result, Report> { match args.alignment_kernel { - AlignmentBackend::Minimap2 => align_with_minimap2_lib(blocks, &args.aln_args), - AlignmentBackend::Mmseqs => align_with_mmseqs(blocks, &args.aln_args), + AlignmentBackend::Minimap2 => align_with_minimap2_lib(blocks, names, &args.aln_args), + AlignmentBackend::Mmseqs => align_with_mmseqs(blocks, names, &args.aln_args), } .wrap_err_with(|| format!("When trying to align sequences using {}", &args.alignment_kernel)) } -pub fn filter_matches(alns: &[Alignment], args: &AlignmentArgs) -> Vec { +pub fn filter_matches(alns: &[Alignment], names: &BlockNames, args: &AlignmentArgs) -> Vec { // - evaluates the energy of the alignments // - keeps only matches with E < 0 // - sorts them by energy @@ -192,11 +199,22 @@ pub fn filter_matches(alns: &[Alignment], args: &AlignmentArgs) -> Vec qry_block.depth() { - AnchorBlock::Ref - } else { - AnchorBlock::Qry - } - } else { - // Equal depth: prefer fewer Ns in the aligned interval (ref wins ties via <=) - let ref_n = ref_block.consensus()[m.reff.interval.to_range()] - .iter() - .filter(|c| c.0 == b'N') - .count(); - let qry_n = qry_block.consensus()[m.qry.interval.to_range()] + + let n_count = |block: &PangraphBlock, interval: &Interval| { + block.consensus()[interval.to_range()] .iter() .filter(|c| c.0 == b'N') - .count(); - if ref_n <= qry_n { - AnchorBlock::Ref - } else { - AnchorBlock::Qry - } + .count() }; - m.anchor_block = Some(anchor); + + let ordering = ref_block + .depth() + .cmp(&qry_block.depth()) + .then_with(|| n_count(qry_block, &m.qry.interval).cmp(&n_count(ref_block, &m.reff.interval))) + .then_with(|| names.sort_key(m.qry.name).cmp(&names.sort_key(m.reff.name))); + + m.anchor_block = Some(match ordering { + Ordering::Less => AnchorBlock::Qry, + _ => AnchorBlock::Ref, + }); } } @@ -408,13 +411,14 @@ fn split_block( pub fn reweave( mergers: &mut [Alignment], mut graph: Pangraph, + names: &BlockNames, thr_len: usize, ) -> Result<(Pangraph, Vec), Report> { // for each merger, assign a new block id (hash of original block ids and alignment intervals) assign_new_block_ids(mergers); - // and decide which block is the anchor, based on depth and n. of ambiguous nucleotides - // (ref block wins ties) - assign_anchor_block(mergers, &graph); + // and decide which block is the anchor, based on depth, n. of ambiguous nucleotides + // and a content-derived tie-break + assign_anchor_block(mergers, &graph, names); // dictionary of BlockId -> alignments. Nb: each alignment is present twice. // this is done to quickly access all alignments for a given block @@ -631,7 +635,8 @@ mod tests { }; let mut mergers = vec![new_aln(1, 2), new_aln(3, 4), new_aln(4, 1)]; - assign_anchor_block(&mut mergers, &pangraph); + let names = BlockNames::from_blocks(&pangraph.blocks); + assign_anchor_block(&mut mergers, &pangraph, &names); assert_eq!(mergers[0].anchor_block, Some(AnchorBlock::Qry)); assert_eq!(mergers[1].anchor_block, Some(AnchorBlock::Ref)); @@ -1005,7 +1010,8 @@ mod tests { let (G, mut M) = generate_example(); let O = G.clone(); let thr_len = 90; - let (G, P) = reweave(&mut M, G, thr_len)?; + let names = BlockNames::from_blocks(&G.blocks); + let (G, P) = reweave(&mut M, G, &names, thr_len)?; // new paths let p1 = &G.paths[&PathId(100)]; @@ -1117,17 +1123,21 @@ mod tests { // CIGAR modified by right 50 bp overhang in ref assert_eq!(p3.cigar, Cigar::from_str("100M50D")?); + // Blocks 30 and 50 have equal depth (2) and no ambiguous bases, so this merger is a full tie. + // The anchor is therefore picked on a content-derived key rather than defaulting to the + // reference, and here it resolves to the query side (block 30). The alignment was recorded as + // qry=30 -> ref=50, so anchoring on the query inverts the cigar: insertions become deletions. let bid50_2 = G.nodes[&nid_100_2].block_id(); let p4 = &p_dict[&bid50_2]; assert_eq!(p4.orientation, Forward); assert_eq!(p4.anchor_block.id(), bid50_2); + assert_eq!(p4.anchor_block.consensus(), &O.blocks[&BlockId(30)].consensus()[0..100]); + assert_eq!(p4.append_block.id(), G.nodes[&nid_300_5].block_id()); assert_eq!( - p4.anchor_block.consensus(), + p4.append_block.consensus(), &O.blocks[&BlockId(50)].consensus()[150..250] ); - assert_eq!(p4.append_block.id(), G.nodes[&nid_300_5].block_id()); - assert_eq!(p4.append_block.consensus(), &O.blocks[&BlockId(30)].consensus()[0..100]); - assert_eq!(p4.cigar, Cigar::from_str("80M10I10M10D")?); + assert_eq!(p4.cigar, Cigar::from_str("80M10D10M10I")?); // assert_eq!(p1.append_block.id(), p1.anchor_block.id()); // assert_eq!(p2.append_block.id(), p2.anchor_block.id()); @@ -1211,8 +1221,9 @@ mod tests { // -- N count tie-breaker (equal depth) -- #[case::equal_depth_ref_fewer_ns (("ATCG", 2), ("NNCG", 2), (2, (0, 4), 1, (0, 4)), AnchorBlock::Ref)] #[case::equal_depth_qry_fewer_ns (("ATCG", 2), ("NNCG", 2), (1, (0, 4), 2, (0, 4)), AnchorBlock::Qry)] - #[case::equal_depth_equal_ns_ref_wins(("ANCG", 2), ("TNCG", 2), (2, (0, 4), 1, (0, 4)), AnchorBlock::Ref)] - #[case::equal_depth_zero_ns_ref_wins(("ATCG", 2), ("GCTA", 2), (2, (0, 4), 1, (0, 4)), AnchorBlock::Ref)] + // Fully-tied cases (equal depth, equal N count) are covered by + // `test_assign_anchor_block_tie_is_invariant_under_role_swap`: they no longer have a fixed + // Ref/Qry answer, since the tie-break is now content-derived rather than reference-first. #[case::equal_depth_many_ns_qry_wins(("NNNG", 2), ("NNCG", 2), (2, (0, 4), 1, (0, 4)), AnchorBlock::Qry)] // -- depth wins over N count -- #[case::qry_deeper_wins (("NNCG", 3), ("ATCG", 2), (1, (0, 4), 2, (0, 4)), AnchorBlock::Qry)] @@ -1268,7 +1279,65 @@ mod tests { align: None, }]; - assign_anchor_block(&mut mergers, &pangraph); + let names = BlockNames::from_blocks(&pangraph.blocks); + assign_anchor_block(&mut mergers, &pangraph, &names); assert_eq!(mergers[0].anchor_block, Some(expected)); } + + /// When depth and ambiguous-base count are both tied, the anchor must be chosen from the block + /// contents, not from which side the aligner happened to call "reference". Swapping the + /// query/reference roles must therefore still select the *same* block. + /// + /// This is the defect that made `pangraph build` sensitive to the order of its input FASTA + /// arguments: `minimap2 -X` picks the direction of each pair by comparing sequence names, which + /// were block ids derived from input position, and the old tie-break then always kept the + /// reference consensus. + #[rstest] + #[case::equal_ns("ANCG", "TNCG")] + #[case::zero_ns("ATCG", "GCTA")] + #[trace] + fn test_assign_anchor_block_tie_is_invariant_under_role_swap(#[case] seq1: &str, #[case] seq2: &str) { + let edits = + |offset: usize| -> BTreeMap { (0..2).map(|i| (NodeId(offset + i), Edit::empty())).collect() }; + + let pangraph = Pangraph { + blocks: btreemap! { + BlockId(1) => PangraphBlock::new(BlockId(1), seq1, edits(0)), + BlockId(2) => PangraphBlock::new(BlockId(2), seq2, edits(100)), + }, + paths: btreemap! {}, + nodes: btreemap! {}, + }; + let names = BlockNames::from_blocks(&pangraph.blocks); + + let hit = |id: usize| Hit { + name: BlockId(id), + length: seq1.len().max(seq2.len()), + interval: Interval::new(0, 4), + }; + let aln = |qry_id: usize, ref_id: usize| Alignment { + qry: hit(qry_id), + reff: hit(ref_id), + matches: 0, + length: 0, + quality: 0, + orientation: Forward, + new_block_id: None, + anchor_block: None, + cigar: Cigar::default(), + divergence: None, + align: None, + }; + + // The same pair of blocks, with the roles swapped. + let mut mergers = vec![aln(2, 1), aln(1, 2)]; + assign_anchor_block(&mut mergers, &pangraph, &names); + + // Resolve each choice back to the block it actually selected. + let anchor_of = |m: &Alignment| match m.anchor_block.unwrap() { + AnchorBlock::Ref => m.reff.name, + AnchorBlock::Qry => m.qry.name, + }; + assert_eq!(anchor_of(&mergers[0]), anchor_of(&mergers[1])); + } } diff --git a/packages/pangraph/src/tree/neighbor_joining.rs b/packages/pangraph/src/tree/neighbor_joining.rs index 8563d54e..c64fe2b7 100644 --- a/packages/pangraph/src/tree/neighbor_joining.rs +++ b/packages/pangraph/src/tree/neighbor_joining.rs @@ -1,5 +1,6 @@ #![allow(non_snake_case)] +use crate::align::block_names::consensus_hash; use crate::distance::mash::mash_distance::mash_distance; use crate::distance::mash::minimizer::MinimizersParams; use crate::pangraph::pangraph::Pangraph; @@ -12,8 +13,21 @@ use itertools::Itertools; use ndarray::{Array1, Array2, Axis, s}; use ndarray_stats::QuantileExt; +/// Content-derived key for ordering the leaves fed to neighbor joining. +/// +/// `Q.argmin()` returns the *first* minimum in row-major order, so the row order of the distance +/// matrix silently decides which pair is joined whenever the Q matrix ties. Keying the order on +/// block content instead of on input position keeps that choice independent of the order the input +/// FASTA files were listed in. +fn graph_content_key(graph: &Pangraph) -> Vec { + graph.blocks.values().map(consensus_hash).sorted().collect() +} + /// Generate guide tree using neighbor joining method. pub fn build_tree_using_neighbor_joining(graphs: Vec) -> Result>>, Report> { + let mut graphs = graphs; + graphs.sort_by_cached_key(graph_content_key); + let mut distances = calculate_distances(&graphs); let mut nodes = graphs diff --git a/packages/pangraph/tests/itest_build_order_invariance.rs b/packages/pangraph/tests/itest_build_order_invariance.rs new file mode 100644 index 00000000..bbeb2e03 --- /dev/null +++ b/packages/pangraph/tests/itest_build_order_invariance.rs @@ -0,0 +1,162 @@ +#[cfg(test)] +mod tests { + use eyre::Report; + use itertools::Itertools; + use pangraph::commands::build::build_args::PangraphBuildArgs; + use pangraph::commands::build::build_run::build; + use pangraph::io::fasta::FastaRecord; + use pangraph::pangraph::pangraph::Pangraph; + use pangraph::representation::seq::Seq; + use pretty_assertions::assert_eq; + use rstest::rstest; + use std::io::Write; + use tempfile::NamedTempFile; + + /// Deterministic pseudo-random nucleotide sequence, so the test does not depend on an RNG. + fn random_seq(len: usize, seed: u64) -> String { + let mut state = seed; + std::iter::repeat_with(|| { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + b"ACGT"[((state >> 33) % 4) as usize] as char + }) + .take(len) + .collect() + } + + /// Applies a substitution every `period` bases, to make sequences homologous but not identical. + fn mutate(seq: &str, period: usize) -> String { + seq + .chars() + .enumerate() + .map(|(i, c)| { + if i % period == 0 { + match c { + 'A' => 'G', + 'G' => 'A', + 'C' => 'T', + _ => 'C', + } + } else { + c + } + }) + .collect() + } + + /// Four genomes sharing two homologous core segments, each with its own accessory segment. + /// + /// Four genomes (rather than three) matter: with a balanced topology the final merge joins two + /// clades of *equal* depth, and with no ambiguous bases the anchor choice is a full tie. That is + /// the case the old reference-first tie-break resolved using the input order, so a three-genome + /// fixture - where depth always decides - cannot detect the defect. + fn genomes() -> Vec<(String, String)> { + let core_a = random_seq(2500, 42); + let core_b = random_seq(2500, 7); + let mosaic = + |a: usize, b: usize, acc: u64| format!("{}{}{}", mutate(&core_a, a), random_seq(700, acc), mutate(&core_b, b)); + vec![ + ("g1".to_owned(), mosaic(50, 70, 1)), + ("g2".to_owned(), mosaic(55, 75, 2)), + ("g3".to_owned(), mosaic(60, 80, 3)), + ("g4".to_owned(), mosaic(65, 85, 4)), + ] + } + + /// Builds a graph from the given genomes, in the given order, optionally with a guide tree. + fn build_in_order(order: &[usize], newick: Option<&str>) -> Result { + let all = genomes(); + let fastas = order + .iter() + .enumerate() + .map(|(index, &g)| FastaRecord { + seq_name: all[g].0.clone(), + desc: None, + seq: Seq::from_str(&all[g].1), + index, + }) + .collect_vec(); + + // The tree file has to outlive the build call. + let tree_file = newick + .map(|nwk| -> Result { + let mut f = NamedTempFile::new()?; + write!(f, "{nwk}")?; + f.flush()?; + Ok(f) + }) + .transpose()?; + + let args = PangraphBuildArgs { + guide_tree: tree_file.as_ref().map(|f| f.path().to_owned()), + ..PangraphBuildArgs::default() + }; + + build(fastas, &args, true) + } + + /// Multiset of block consensus sequences: the graph structure, independent of block identifiers. + fn consensus_multiset(graph: &Pangraph) -> Vec { + graph + .blocks + .values() + .map(|b| b.consensus().as_str().to_owned()) + .sorted() + .collect() + } + + /// All six argument orders must produce the same graph. + /// + /// Block identifiers are derived from the position of a record in the input, and used to be + /// passed to the aligner as sequence names. Since `minimap2 -X` picks the direction of each + /// pairwise alignment by comparing those names, argument order silently decided which block was + /// query and which was reference - and the reference's consensus was the one the merged block + /// inherited. + /// + /// The `no_guide_tree` case additionally guards the neighbor-joining path: `Q.argmin()` resolves + /// ties by matrix row order, so unless the leaves are ordered by content the inferred topology + /// itself flips with the argument order. + #[rstest] + #[case::no_guide_tree(None)] + #[case::balanced_guide_tree(Some("((g1,g2),(g3,g4));"))] + #[case::ladder_guide_tree(Some("(((g1,g2),g3),g4);"))] + #[trace] + fn test_build_is_invariant_under_argument_order(#[case] newick: Option<&str>) -> Result<(), Report> { + let permutations = [ + [0, 1, 2, 3], + [1, 0, 2, 3], + [3, 2, 1, 0], + [2, 3, 0, 1], + [1, 3, 0, 2], + [3, 0, 2, 1], + ]; + + let expected = consensus_multiset(&build_in_order(&permutations[0], newick)?); + assert!( + expected.len() > 1, + "expected a non-trivial graph, got {} block(s)", + expected.len() + ); + + for order in &permutations[1..] { + let actual = consensus_multiset(&build_in_order(order, newick)?); + assert_eq!(expected, actual, "argument order {order:?} produced a different graph"); + } + + Ok(()) + } + + /// Swapping two siblings in the guide tree does not change the topology, so it must not change + /// the graph either. `graph_join` is symmetric and `merge_graphs` uses left/right only for + /// logging, so this holds once the aligner input is canonical. + #[rstest] + fn test_build_is_invariant_under_guide_tree_sibling_swap() -> Result<(), Report> { + let straight = build_in_order(&[0, 1, 2, 3], Some("((g1,g2),(g3,g4));"))?; + let swapped = build_in_order(&[0, 1, 2, 3], Some("((g2,g1),(g4,g3));"))?; + + assert_eq!(consensus_multiset(&straight), consensus_multiset(&swapped)); + + Ok(()) + } +}