diff --git a/.gitea/workflows/ci.yml b/.gitea/workflows/ci.yml index 1f8e17d..dfaf01e 100644 --- a/.gitea/workflows/ci.yml +++ b/.gitea/workflows/ci.yml @@ -25,8 +25,14 @@ jobs: ~/.cargo/registry ~/.cargo/git src/target - key: ${{ runner.os }}-cargo-${{ hashFiles('src/Cargo.lock') }} - restore-keys: ${{ runner.os }}-cargo- + # v2: bump this salt to force a clean cache when a stuck/killed job + # may have saved a corrupted incremental-compilation `src/target` + # (observed 2026-08-11: a stale test binary deadlocked in NUMA + # worker startup; a from-scratch rebuild of the same source fixed + # it instantly, pointing at cache corruption rather than a source + # bug). + key: ${{ runner.os }}-cargo-v2-${{ hashFiles('src/Cargo.lock') }} + restore-keys: ${{ runner.os }}-cargo-v2- - name: Build run: cargo build --release diff --git a/docmd/theory/evolutionary_distances.md b/docmd/theory/evolutionary_distances.md index ba63278..8c0e833 100644 --- a/docmd/theory/evolutionary_distances.md +++ b/docmd/theory/evolutionary_distances.md @@ -203,6 +203,180 @@ implication is that columns should be allowed partial coverage (>=2 resolved genomes, not unanimous) rather than requiring every genome to be net single-copy at that locus. +## Context, detectability, and a 3-way ordinal distance per pair + +Empirical follow-up to the pseudo-alignment idea above: `obikmer distance +--snp` was run on a real 20-genome benchmark index and the resulting FASTA +fed to `raxml-ng`. Two problems surfaced, both traced back to conflating +distinct notions under one symbol. + +**"Context", precisely.** Sharing a central base between two genomes is not +just sharing a nucleotide — it is sharing a **context**: the `2m` flanking +bases, identical, which is a homology claim about that flanked window +(guaranteed non-coincidental by k-specificity, Bias 4 above), *not* a claim +about orthology or paralogy of the copy each genome carries. `A` opposite `C` += same context, divergent centre. `A` opposite nothing = **this context is +not observed in one of the two genomes** — informative, not neutral. + +**Why the IUPAC/DNA encoding used for the first `--snp` test was wrong.** +Feeding IUPAC-coded ambiguity into a standard DNA model (`raxml-ng --model +GTR+G`) is a semantic mismatch: Felsenstein-pruning ML treats an ambiguous +tip as "exactly one true state, unknown which" (a uniform partial-likelihood +vector over compatible bases), not "these states are simultaneously +present". The two encodings look identical (same IUPAC letters) but the +software reads them backwards from what was intended — this invalidates the +literal branch lengths from that first experiment (topology-level groupings +by genus were still informative, see the worked example further down). + +**Why `-` (absence) must not be scored as similarity, but also must not be +scored as a shared character between two absences.** Two genomes both +lacking a context are not observed to resemble each other at that locus — +neither is observed to differ from the other either. It is a symmetric +non-observation, uninformative for that pair, and should contribute nothing +(not a small positive nor a small negative signal) to their distance. A +genome carrying a state (`A`) against one carrying none is a different case +entirely: informative, and should not be scored as neutral "missing data" +the way a generic DNA/ML pipeline would. + +**Detectability vs existence — a deliberate simplification, accepted.** +"Context not observed" conflates two different biological events: (1) true +loss of the locus, (2) the locus still exists but a mutation/indel *outside* +the centre, anywhere in the `2m` flanks, broke k-mer recognition. The design +adopts a rigorist stance on purpose: any flank-breaking mutation counts as +"this context no longer exists", full stop — because both causes (1) and (2) +independently require *at least* one more mutational event than a lone +central substitution would. This licenses treating "context absent in one of +the two genomes" as a **lower bound** on distance strictly greater than a +plain central SNP, without needing to know which of the two causes applies. +Coarser than a true event count, and accepted as such (fine substitution-type +modelling, e.g. transition/transversion weighting, is a secondary +refinement, not required for this to be useful). + +**Resulting ordinal distance between two genomes at one family/context:** + +| Comparison | Distance | Meaning | +|---|---|---| +| same centre (`A`/`A`) | `0` | identical | +| different centre, both single-copy (`A`/`C`) | `1` | plain central SNP | +| one genome has a state, the other has none | `>1` (lower bound) | context undetectable in one genome — at least one extra mutational event, of unknown type | +| neither genome has any state (`∅`/`∅`) | excluded | symmetric non-observation, not comparable, contributes nothing | + +This is a direct extension of `KmerIndex::raw_snp_distance` (`obikindex/src/siblings.rs`), +which today only implements the `0`/`1` rows and silently drops everything +else (including the informative `>1` row) rather than scoring it. + +**Open, not yet resolved:** +- Calibrating `>1` to a real number for tools expecting continuous distances + (NJ/UPGMA, ML branch lengths), rather than an arbitrary placeholder. + Natural route: estimate `p_hat` from the resolved (`0`/`1`) sites first, + then use the already-derived ascertainment formula (`P(usable window + showing a central SNP) = p * (1-p)^(2m)`, Bias 1 above) to derive a + model-consistent value for the `>1` bucket instead of guessing a constant. +- Where multi-copy/ambiguous states (the IUPAC case: a genome carrying more + than one form) fit into this ordinal scheme — plausibly also `>1` by the + same "at least one extra event" argument (a second form appearing is a + gain, itself an event), but not yet worked out. + +**Practical alternative validated for the pseudo-alignment output itself** +(orthogonal to the ordinal-distance question above, useful regardless of how +`>1` ends up calibrated): re-encode each family as 4 independent binary +presence/absence characters (`A`,`C`,`G`,`T` columns) instead of one IUPAC +column, feed to a `BIN`-type model instead of `DNA`. `∅` becomes an explicit +`0000` state (identity with another `0000`, not missing data) rather than a +gap — removes the semantic mismatch above by construction. Known cost, +accepted for now: a plain central substitution (`A` -> `C`) becomes 2 binary +flips (`1000` -> `0100`), overweighting substitutions relative to true +gain/loss events, and the 4 sub-characters of one family are not +statistically independent the way a generic `BIN` model assumes. A proper +fix (single 16-state alphabet, i.e. the powerset of `{A,C,G,T}`, with a +substitution-rate structure that respects the subset lattice rather than a +fully general 16x16 GTR-analogue) is very likely not expressible in +`raxml-ng`'s `MULTI` datatype as-is (Mk or fully-general rates only) and a +fully general 16-state rate matrix is almost certainly unidentifiable here +(states of cardinality >=3 are ~2% of sites in the benchmark run). Treated as +a longer-term research question, not a near-term implementation target. + +## Sankoff parsimony as the resolution of the 16-state model problem + +The "longer-term research question" just above (a 16-state alphabet — the +powerset of `{A,C,G,T}` — with a substitution structure that respects the +subset lattice) turns out to have a near-term answer, once the *unification* +question below is worked through. + +**Distance methods (NJ/UPGMA/ME) vs. character methods (parsimony/ML): not a +deep philosophical divide, but a real practical distinction for this +project.** Historically "phenetic" (characters -> distances -> tree) and +"cladistic" (characters -> tree directly) approaches were presented as +opposed schools; the modern view is a mathematical continuity, not a +dichotomy — Minimum Evolution (ME: find the tree minimizing total branch +length from a distance matrix) and Maximum Parsimony (MP: find the tree +minimizing total character-state changes) are both instances of "minimise a +global explanatory cost", and coincide under simple encodings (see Farris +1983, "The Logical Basis of Phylogenetic Analysis"; the MP/ME connection is +developed in the Minimum Evolution / Balanced Minimum Evolution literature, +e.g. Nei and colleagues — citations not independently re-verified here, flag +before quoting further). NJ's own agglomeration step already uses the whole +distance matrix jointly (the Q-matrix), not just the pair being merged — an +earlier claim in this discussion that distance methods are "blind" to +cross-taxon structure at every stage was too strong. + +What *does* remain a real, structural distinction for this project: in a +character method, a given character's cost is **re-evaluated per candidate +topology** during tree search (the same family can cost 1 change under one +topology, 2 under another). In a pairwise-distance pipeline +(`raw_snp_distance` as it exists today), each family's contribution to +`d(i,j)` is computed **once**, independent of any candidate topology, before +NJ/UPGMA ever runs — so a question like "does this shared `∅` look like a +synapomorphy under topology T" can never be posed in that pipeline, for any +T. That question is only answerable by a method that tests candidate +topologies and re-scores characters under each — i.e. a character method. + +**Sankoff parsimony directly resolves the `∅`/gain-loss/substitution +question, without the identifiability problem of a fitted 16-state model.** +Sankoff's algorithm generalises Fitch parsimony to an arbitrary +user-supplied cost matrix between states (`obikseq`/`obikindex` would treat +each family as a `2^{4}`-state character, state = subset of `{A,C,G,T}` +observed in that genome, `∅` included as a real state, not a gap). The +previous 16-state idea failed specifically because *fitting* a full 16x16 +rate matrix by ML is unidentifiable at this data volume; Sankoff sidesteps +that because the cost matrix is **fixed a priori from domain knowledge**, not +estimated — e.g. `c({A},{C}) = 1` (a substitution), `c({A},{A,C}) = 1` (a +gain), `c({A,C},{A}) = 1` (a loss), `c({A,C},{G,T}) = 2` (two changes) — no +estimation, no overparameterisation. This reframes gain/loss and central +substitution as two cost categories with independently chosen weights, +exactly the "two families of parameters" (`mu_substitution`, `mu_gain/loss`) +floated earlier in this discussion, now with an actual algorithmic home. + +**Caveat, not blocking for this project's scope.** Sankoff is still +parsimony: in principle exposed to Felsenstein's statistical-inconsistency +result under long-branch attraction (already invoked earlier against +`D_F = min(a,b)`) — parsimony and ML only provably coincide in the +short-branch regime. This is not a practical concern here because it is +exactly this estimator's declared target (closely related genomes, short +branches) — the regime where parsimony's known failure mode does not apply — +but worth stating explicitly as a scope guard rather than leaving it +implicit. + +**Cheapest next experiment: don't write a Sankoff tree-search engine, use +one that exists.** The hard part of a from-scratch implementation is not the +Sankoff DP itself (a straightforward dynamic program over a *fixed* tree) but +the topology search (SPR/NNI with incremental re-scoring) that comes for +free with `raxml-ng` on the ML side. **TNT** (Tree analysis using New +Technology, free, standard in morphological cladistics) already implements +Sankoff parsimony with a custom cost matrix plus topology search — the +family-state matrix (already close to what `--snp` produces, minus the +IUPAC/DNA-model mismatch) could be fed there directly, no new code required, +before considering a bespoke engine. + +**The concrete comparison this unlocks:** run both pipelines on the same +family data — +`k-mer families -> D_ij -> NJ/ME/UPGMA` (phenetic, what exists today) vs. +`k-mer families -> characters -> argmin_T Sankoff-cost(T)` (cladistic, via +TNT) — and compare the resulting topologies. Agreement would validate that +the pairwise-distance projection preserves the phylogenetic signal; +disagreement would pinpoint exactly what the projection to a single number +per pair loses. Not yet run. + ## Heterozygosity, ploidy, and consensus-assembly inputs A within-genome multiplicity signal (more than one of the 4 central forms diff --git a/src/Cargo.lock b/src/Cargo.lock index a4703b9..2318aa7 100644 --- a/src/Cargo.lock +++ b/src/Cargo.lock @@ -1715,7 +1715,7 @@ dependencies = [ [[package]] name = "obikmer" -version = "1.1.41" +version = "1.1.42" dependencies = [ "clap", "csv", diff --git a/src/obikmer/Cargo.toml b/src/obikmer/Cargo.toml index 49e2485..d86ab0e 100644 --- a/src/obikmer/Cargo.toml +++ b/src/obikmer/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "obikmer" -version = "1.1.41" +version = "1.1.42" edition = "2024" [[bin]]