From c26623fa00a990884883d3ecd5ddd374cbc2f27d Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Wed, 12 Aug 2026 19:51:35 +0200 Subject: [PATCH] Add repeatable --exclude-genome flag to obikmer distance command Integrate in-memory row/column zeroing and alignment filtering across SNP, Sankoff, TNT, and IQ-TREE output paths. Add strict validation for missing labels, implement `write_iqtree` with empirical stationary frequencies, and introduce `--raw-snp-counts` diagnostic CSV output. Update documentation to reflect experimental validation of backbone resolution limits and theoretical considerations for CTMC rate matrices. --- docmd/theory/evolutionary_distances.md | 207 ++++++++++++++++++++++++- src/obikmer/src/cmd/distance/mod.rs | 83 +++++++++- 2 files changed, 282 insertions(+), 8 deletions(-) diff --git a/docmd/theory/evolutionary_distances.md b/docmd/theory/evolutionary_distances.md index 6fd21d8b..c7e2ebfe 100644 --- a/docmd/theory/evolutionary_distances.md +++ b/docmd/theory/evolutionary_distances.md @@ -897,7 +897,212 @@ Models](https://iqtree.github.io/doc/Complex-Models). (see "A concrete Sankoff cost matrix" above) is already a log-rate — a genuine CTMC rate matrix `Q` could plausibly be recovered as `rate(a,b) = exp(-cost(a,b))`, renormalised so each row sums to zero, once -the stationary-frequency gap above is closed. Not yet attempted. +the stationary-frequency gap above is closed. Superseded by the two +sections below, which implement and then substantially revise this. + +### `R`/`π` implemented; `--exclude-genome` added; rogue-taxon test negative (2026-08-12) + +`write_iqtree` (`obikmer/src/cmd/distance/iqtree.rs`) implements exactly +the `R = exp(-cost)` / empirical-`π` design above: writes +`_iqtree.model` (lower-triangular `R`, PAML order, then `π`) and +`_iqtree.fasta` (alignment recoded to the compact `0..k-1` +alphabet of states actually present), states/frequencies restricted to +whichever of the 16 canonical states actually occur, compactly +renumbered — see the `MORPH{N}` risk above. Verified end to end on the +20-genome benchmark: real, non-zero, varied branch lengths (`Total tree +length: 6.773`), converged log-likelihood, `State frequencies:` in +IQ-TREE's own output matching the computed `π` exactly. + +**Rogue-taxon hypothesis, tested and refuted.** The backbone (*Yersinia*, +*Proteus*, *Opitutus*, *Shouchella*, *Wolbachia*...) resolves as a +near-linear comb with ~8 near-zero branch lengths — visible on the real +tree, and independently on a Mash+NJ tree built from an entirely +different signal (whole-genome k-mer distance, no relation to the Sankoff +pipeline), which shows the same comb shape. Hypothesis: `Saccharolobus`/ +`Candidozyma` (near-zero real signal — see below) destabilise the +heuristic tree search enough to also blur resolution elsewhere ("rogue +taxa", a documented phenomenon). Tested directly: reran `iqtree3` after +removing both taxa (`--exclude-genome`, `+ASC` recomputed on the +resulting variable-sites-only alignment since removing taxa turns some +columns invariant) — **still exactly 8 near-zero backbone branches**, +identical comb shape. Refuted for this dataset: the backbone's weak +resolution is a property of the character system's signal at that +divergence depth (matches "Run 3" above), not rogue-taxon interference. + +**Diagnosed why `Saccharolobus`/`Candidozyma` place so poorly**, using a +new diagnostic (`--raw-snp-counts`, `_rawsnp_counts.csv`: one row +per genome pair, `n_snp,n_shared,n_eligible,ratio` — the counts +`--raw-snp-distance`'s ratio-only matrix discards, needed because +`ratio=0.0` from 2 eligible loci and from 2000 look identical in the +ratio alone). `Saccharolobus` has 21,020 non-`∅` sites in the real +alignment (comparable to other taxa) — but **21,017 of them (100.0%) are +private**: no other genome has a non-`∅` state at the same site. Only 3 +sites are shared with any other genome at all (1, 2, and 14 +co-occurring genomes respectively). Real informativeness for placement +tracks shared sites, not raw non-`∅` count — with no other archaeon (or +even archaea-adjacent bacterium) in the dataset, there is nothing to +anchor `Saccharolobus`'s position to, regardless of how much of its own +data exists. `Candidozyma` shows the same pattern, more extreme (almost +all `NA`, the 3 non-`NA` pairs all exactly `0.0`, never `1.0` — itself a +tell: with `n_eligible=1`, the ratio can only be exactly `0` or `1`, so +3-for-3 landing on `0.0` is more than sampling noise alone would predict; +possibly ascertainment bias — see below — or possibly a few genuinely +ultra-conserved loci; not resolved). + +**`--exclude-genome LABEL`** (repeatable, `obikmer distance`) added for +exactly this kind of test: zeroes the excluded genome's row/column in +`RawSnpDistanceOutput` after `raw_snp_distance` runs (a pair with zero +counts is already skipped by `calibrate_p_hat`/`base_pair_tally` — no +`obikindex` traversal change needed) and drops its row from +`SnpAlignment` before any output is written. Deliberately *not* index +surgery (a new, smaller on-disk index) — genome sets to exclude are +expected to change between quick tests, so an in-memory filter is the +right tool, not a new index-rewriting subsystem. Scoped to the sibling-annex +family of computations (`--raw-snp-distance`/`--raw-snp-counts`/`--snp`/ +`--sankoff` and everything it implies) — does not affect the plain +`--metric` distance matrix/NJ/UPGMA path (a different, unrelated +computation on `idx.distance()`, not touched). + +**Caveat surfaced while reusing a stale `π`**: rerunning IQ-TREE on a +genome-reduced alignment while keeping the *original* (20-genome) +model file is inconsistent — `π`'s composition shifts once low-cardinality-0 +columns that were only variable because of the removed taxa drop out. +Measured directly: `π(∅)` `0.9048` (20 genomes) → `0.8869` (18 genomes, +variable sites only) — real (~378k affected cells, from `Saccharolobus`'s +~21k private sites × 18 remaining genomes) but modest in *proportion* +(~2 points) because the alignment was already monomorphic-filtered +before `Saccharolobus` was ever added, so removing it only drops the +subset of columns that were variable *because of* it specifically, not +every column it appears in. + +### `R` via `exp(-cost)` is wrong for a CTMC; cardinality/composition decomposition (open, 2026-08-12) + +**The flaw in `R = exp(-cost)`, precisely.** `cost` (`build_cost_matrix`'s +output) is a *shortest-path closure* over an elementary-edit graph +(Floyd-Warshall) — correct and required for Sankoff parsimony, where +`cost(a,b)` must be a metric. But a CTMC's own matrix exponential +(`exp(Qt)`, computed internally by IQ-TREE) *already* sums over +paths of every length through the elementary rates — that's the whole +mechanism by which a CTMC generates indirect transitions. Feeding it a +pre-summed, multi-hop shortest-path cost and exponentiating that *again* +as if each entry were a direct edge double-applies the "compose multiple +steps" logic once in log-space (Floyd-Warshall, additive) and once more +inside IQ-TREE's own exponential — systematically over-penalising +non-adjacent state pairs (e.g. `∅→{A,C}` priced as two chained edges, +`2×c_ctx`, when it should be one direct lookup). + +**Resolution, in two parts — both estimated directly from the real +alignment, not smoothed through a small parametric formula:** + +1. **Cardinality model**: a 5-state (`0,1,2,3,4`) first-order Markov + chain, estimated from the empirical cardinality co-occurrence table + (pooled across all *included* genome pairs — same + saturated/`NA`-pair exclusion discipline as `calibrate_p_hat`, 14 + saturated + 18 `NA` pairs excluded of 190 in the benchmark; the + result is materially different from the unfiltered version and more + internally consistent, not just "cleaner"). Diagonal included (the + probability of a family *staying* at the same cardinality is part of + the model, not assumed away). Measured, cost `= -ln(observed/expected + under independence)`, on the 20-genome benchmark (`π` here from the + marginal cardinality distribution: `90.481%, 5.261%, 4.183%, 0.067%, + 0.008%` for `c=0..4` respectively — note `c=2` is *not* rare, almost + as common as `c=1`): + + ``` + c=0/c=0: 0.166 c=0/c=1: 0.160 c=0/c=2: 0.157 c=0/c=3: 0.059 c=0/c=4: 0.006 + c=1/c=1: -0.803 (enriched) c=1/c=2: 0.381 c=1/c=3: 1.027 c=1/c=4: 1.794 + c=2/c=2: 4.032 (sharply suppressed) c=2/c=3: 2.620 c=2/c=4: 2.416 + c=3/c=3: -0.254 c=3/c=4: -0.411 (too few observations to trust) + ``` + Two robust findings, both stable under the saturation filter: (a) all + `∅`-involving costs are low and close to each other (`0.006`–`0.166`) + regardless of how many members are gained/lost at once — sharply at + odds with the current model's implicit `×2`/`×3` multi-hop scaling; + (b) `c=2/c=2` (two genomes both showing an ambiguous 2-member state at + the same site) is dramatically under-represented (~1.5–1.8% of the + independence expectation) — a real, robust anomaly, not explained. + +2. **Composition model**: unchanged — `sub_cost` (the existing 6-category + Ts/Tv-biased calibration), estimated **only from unambiguous sites** + (cardinality-1 ↔ cardinality-1 pairs), because composition bias is a + substitution phenomenon and only means something when cardinality is + conserved. Checked whether composition bias also appears in pure + gain/loss events (no substitution involved, so no bias expected a + priori): single-base "gain" (`∅→`single base) is close to uniform + (`24.35/25.92/26.22/23.51%` for A/C/G/T) — consistent with "no + mutational mechanism, so no bias" as expected. Two-base "gain" + (`∅→`2-member state) is *not* uniform even after correcting for the + real (non-25/25/25/25) marginal base frequencies: `{A,G}`/`{C,T}` + (the transition-linked pairs) mildly enriched (`obs/exp` `1.12`/ + `1.16`), `{C,G}` sharply suppressed (`obs/exp 0.563`, expected to be + the *most* common pair under independence since C and G are + individually the two most frequent bases, observed the least) — an + unexplained anomaly, deliberately **not** built into the model (no + mechanistic story for why gain/loss would carry a `{C,G}`-specific + bias), left as an open puzzle rather than fit. + +**Composing the two into a full pairwise cost — the actual replacement +for both the elementary-edge graph and its Floyd-Warshall closure.** For +any two states `A`, `B` (not just the "clean" same-cardinality or +pure-subset cases the current graph handles directly): let `shared = A∩B` +(free), `lost = A\B`, `gained = B\A`. Pair off `k = min(|lost|,|gained|)` +elements between `lost` and `gained` as substitution events, choosing the +pairing that minimises total `sub_cost` (a trivial assignment problem — +at most 4 elements per side). The `|lost|-k` (or `|gained|-k`) leftover, +unpaired elements are a *pure* cardinality change, `|A|→|B|`, costed by +**one direct lookup** in the cardinality model above — no chaining, no +Floyd-Warshall. Example: `{A,C}→{G}` (cardinality 2→1, no shared base): +pair 1 substitution (cheaper of `A→G`, `C→G`), 1 base left over unpaired +→ cost `= sub_cost(chosen pair) + cardinality_cost(2→1)`. This is a +direct, closed-form cost for *every* pair of the 16 states, replacing +`build_cost_matrix`'s graph-plus-shortest-path construction outright — +and specifically fixes the CTMC double-counting problem, since every +entry is now a single decomposed lookup, never a sum of chained edges. + +**Reframed as a likelihood (product of probabilities), not a cost (sum of +`-ln`s) — same content, but forces the diagonals (no-change cases) to be +kept rather than implicitly dropped.** Both sub-models are proper +transition *probability* matrices, diagonal included — `P_cardinality` +includes "stay at the same cardinality", `P_composition` includes "stay +the same base" (e.g. `P(G→G)`, not just `sub_cost`'s off-diagonal +entries). Composing: + +``` +P(A→B) = P_cardinality(|A|→|B|) + × ∏_{x ∈ A∩B} P_composition(x→x) (shared bases: "stayed") + × ∏_{(x,y) paired} P_composition(x→y) (parsimony-paired substitutions) +``` + +with the unpaired leftover `lost`/`gained` elements (if `|lost|≠|gained|`) +contributing *nothing further* beyond the `P_cardinality` term already +counted — consistent with the finding above that pure gain/loss carries +no separate composition bias worth modelling. Worked example, +`{A,C}→{A,G}` (`shared={A}`, one paired substitution `C→G`, nothing left +over): `P = P_cardinality(2→2) · P_composition(A→A) · P_composition(C→G)`. + +The row-wise product of these two independently-calibrated models isn't +guaranteed to already sum to exactly `1` across all `B` for a fixed `A` +(the two aren't perfectly independent in reality) — so each row of the +resulting 16×16 matrix is renormalised (divided by its own sum) after +composition, **not** the matrix as a whole (which would produce a joint +distribution over `(A,B)` pairs, the wrong object — a transition matrix +needs each row, "given I start in `A`", to be a valid distribution over +where I end up). + +This now covers every pair of the 16 states with no unhandled case +identified. Gives 120 parameters, but derived from two small, +well-estimated pieces (a 5×5 cardinality model, a 4×4 composition model) +rather than fit or smoothed independently per pair. + +**Status: designed, not implemented.** Would replace `build_cost_matrix` +(`obikindex/src/sankoff.rs`) and the single `c_ctx` scalar/parameter +entirely; `sub_cost`'s own calibration is untouched. Not yet decided +whether the "clean" parsimony-graph cost (`build_cost_matrix`'s current +output, still needed for `--tnt`/`--phyg`) should be replaced too, kept +as a separate simpler approximation, or derived as a special case of this +same decomposition (it likely can be — the same `shared`/`lost`/`gained` +logic, with the cardinality model reduced back to a single scalar `c_ctx` +substituted in, recovers exactly the current graph). ## Heterozygosity, ploidy, and consensus-assembly inputs diff --git a/src/obikmer/src/cmd/distance/mod.rs b/src/obikmer/src/cmd/distance/mod.rs index 3dfe08d3..c3b1035e 100644 --- a/src/obikmer/src/cmd/distance/mod.rs +++ b/src/obikmer/src/cmd/distance/mod.rs @@ -85,6 +85,24 @@ pub struct DistanceArgs { #[arg(long)] pub sibling_annex: bool, + /// Exclude a genome (by its exact label) from every computation below + /// that reads the sibling annex — `--raw-snp-distance`/`--raw-snp-counts`, + /// `--snp`, and `--sankoff` (and everything `--sankoff` implies: + /// `p_hat`, `sub_cost`, `c_ctx`, the exported matrix/alignment, + /// `--tnt`/`--phyg`/`--iqtree`). Repeatable. Does *not* affect the plain + /// `--metric` distance matrix/NJ/UPGMA path (a different, unrelated + /// computation). Applied by zeroing the excluded genome's row/column + /// after `raw_snp_distance` runs (a pair with zero counts is already + /// skipped by `calibrate_p_hat`/`base_pair_tally`, so this needs no + /// change to the underlying traversal) and by dropping its row from + /// `snp_pseudo_alignment`'s output — the annex is still built/scanned + /// for the excluded genome too, just not used afterward. For a genome + /// with almost no informative sites shared with anything else (see + /// `docmd/theory/evolutionary_distances.md`, the IQ-TREE/Mash rogue-taxon + /// discussion), its presence can otherwise silently bias `p_hat`/`R`. + #[arg(long = "exclude-genome", value_name = "LABEL")] + pub exclude_genome: Vec, + /// Tally the sibling-count distribution (CSV) of an already-built annex /// (run with `--sibling-annex` first, in this invocation or an earlier /// one). A separate, occasional diagnostic pass — not run every time the @@ -211,6 +229,52 @@ pub fn run(args: DistanceArgs) { let labels: Vec = idx.meta().genomes.iter().map(|g| g.label.clone()).collect(); let n = labels.len(); + // ── Genome exclusion (`--exclude-genome`) ─────────────────────────────── + // Applied by zeroing a `RawSnpDistanceOutput`'s excluded rows/columns + // (`zero_excluded_pairs`) — `calibrate_p_hat`/`base_pair_tally` already + // skip any pair with zero total counts, so this needs no change to + // `obikindex`'s traversal — and by dropping the excluded genome's row + // from a `SnpAlignment` plus the matching label (`drop_excluded`), + // since an all-`∅` row for an "excluded" genome would otherwise still + // reach TNT/PhyG/IQ-TREE as a real (empty) taxon. + let exclude_mask: Vec = { + let mut mask = vec![false; n]; + for label in &args.exclude_genome { + match labels.iter().position(|l| l == label) { + Some(i) => mask[i] = true, + None => { + eprintln!("error: --exclude-genome {label:?} does not match any genome in this index"); + std::process::exit(1); + } + } + } + mask + }; + let zero_excluded_pairs = |result: &mut RawSnpDistanceOutput| { + for i in 0..n { + if !exclude_mask[i] { + continue; + } + for j in 0..n { + result.snp[[i, j]] = 0; + result.snp[[j, i]] = 0; + result.shared[[i, j]] = 0; + result.shared[[j, i]] = 0; + } + } + }; + let drop_excluded = |alignment: SnpAlignment| -> (SnpAlignment, Vec) { + let sequences = alignment.sequences.into_iter().enumerate() + .filter(|(i, _)| !exclude_mask[*i]) + .map(|(_, seq)| seq) + .collect(); + let kept_labels = labels.iter().enumerate() + .filter(|(i, _)| !exclude_mask[*i]) + .map(|(_, l)| l.clone()) + .collect(); + (SnpAlignment { sequences }, kept_labels) + }; + // ── Sibling-count/minorant annex (independent of the distance metric) ── // Construction (`--sibling-annex`) and stats (`--sibling-stats`) are // deliberately decoupled: the annex is meant to be (re)built routinely, @@ -237,17 +301,19 @@ pub fn run(args: DistanceArgs) { write_sibling_stats_csv(&stats, &labels, &args.output); } if args.raw_snp_distance { - let result = idx.raw_snp_distance().unwrap_or_else(|e| { + let mut result = idx.raw_snp_distance().unwrap_or_else(|e| { eprintln!("error computing raw SNP distance: {e}"); std::process::exit(1); }); + zero_excluded_pairs(&mut result); write_raw_snp_distance_csv(&result, &labels, &args.output); } if args.raw_snp_counts { - let result = idx.raw_snp_distance().unwrap_or_else(|e| { + let mut result = idx.raw_snp_distance().unwrap_or_else(|e| { eprintln!("error computing raw SNP distance: {e}"); std::process::exit(1); }); + zero_excluded_pairs(&mut result); write_raw_snp_counts_csv(&result, &labels, &args.output); } if args.snp { @@ -255,13 +321,15 @@ pub fn run(args: DistanceArgs) { eprintln!("error computing SNP pseudo-alignment: {e}"); std::process::exit(1); }); - write_snp_fasta(&alignment, &labels, &args.output); + let (alignment, kept_labels) = drop_excluded(alignment); + write_snp_fasta(&alignment, &kept_labels, &args.output); } if args.sankoff || args.tnt || args.phyg || args.iqtree { - let raw = idx.raw_snp_distance().unwrap_or_else(|e| { + let mut raw = idx.raw_snp_distance().unwrap_or_else(|e| { eprintln!("error computing raw SNP distance: {e}"); std::process::exit(1); }); + zero_excluded_pairs(&mut raw); let estimate = calibrate_p_hat(&raw, args.sankoff_ratio_ceiling); let m = (idx.kmer_size() - 1) / 2; @@ -282,16 +350,17 @@ pub fn run(args: DistanceArgs) { eprintln!("error computing SNP pseudo-alignment: {e}"); std::process::exit(1); }); - write_sankoff_alignment_fasta(&alignment, &labels, &args.output); + let (alignment, kept_labels) = drop_excluded(alignment); + write_sankoff_alignment_fasta(&alignment, &kept_labels, &args.output); if args.tnt { - write_sankoff_tnt(&matrix, &alignment, &labels, &args.output, args.sankoff_cost_scale); + write_sankoff_tnt(&matrix, &alignment, &kept_labels, &args.output, args.sankoff_cost_scale); } if args.phyg { write_sankoff_phyg(&matrix, &args.output, args.sankoff_cost_scale); } if args.iqtree { - write_iqtree(&matrix, &alignment, &labels, &args.output); + write_iqtree(&matrix, &alignment, &kept_labels, &args.output); } }