obikmer/DevDocMD/theory/evolutionary_distances.md

Central-position SNP distance (discussion)

Not implemented. Design discussion for a substitution-rate estimator that observes SNPs directly from paired-genome k-mer comparison, as an alternative to Mash's Poisson-Jaccard inference (see obicompactvec for the implemented Jaccard/Mash distances).

Motivation

Primary intent: restrict the comparison to what is actually comparable. Mash's Jaccard is computed over the union of both genomes' k-mer content: anything not identically shared is folded into a single undifferentiated mass, whether the cause is a point substitution, a genuinely absent homologous region (lineage-specific content, gene-family expansion, HGT, genome-size asymmetry), or a diverged paralogous copy. The model then back-infers a single mutation rate from that mass, silently attributing non-homology to mutation. The central-SNP approach instead conditions every comparison on local, positive evidence of homology: a locus only enters the statistic if its 2m flanking bases (m = (k-1)/2) are found intact in both genomes — genuinely absent or non-homologous content is excluded from the comparison entirely (neither numerator nor denominator), rather than silently counted as divergence. This is a conditioning on comparability, not just a richer summary statistic — see "Statistic and correspondence with shared" below for how it plays out against genome-size asymmetry and diverged gene families, and "Heterozygosity, ploidy, and consensus-assembly inputs" for the corresponding paralogy/heterozygosity filter.

Secondary benefit: access to the substitution's nature. Because the central base of an odd-k window is directly observable once the flanks are confirmed conserved, this also yields more than a rate — the transition/transversion split — enabling classical corrected distances (Jukes-Cantor, Kimura 2-parameter, LogDet) that a single Jaccard scalar cannot support.

Statistic and correspondence with shared

A genomic position p is covered by k overlapping k-mer windows. Requiring the substitution to sit at the window's center makes exactly one window per SNP eligible — a 1:1 correspondence between SNP and center-neighbor k-mer pair, avoiding the ~k-fold overcount of an any-position neighbor search.

A locus with a fully conserved k-window (flanks and center) is an exact-shared k-mer at that locus; a locus with conserved flanks but a substituted center is a "central SNP". Both count each locus exactly once, in matching units:

p_hat[i,j] = SNP[i,j] / (SNP[i,j] + shared[i,j])

p_hat is P(center substituted | 2m flanks conserved). shared[i,j] here is not the general-purpose shared_kmers matrix used by Jaccard/Mash (--shared-kmers, BitPartials::partial_jaccard / CountPartials::partial_threshold_jaccard) — that matrix counts raw k-mer identity with no per-genome copy-number constraint, whereas p_hat's denominator applies the eligibility rule defined below (raw or paralogy-filtered). Both SNP and shared are accumulated by the same sweep, from the same per-locus candidate set (source k-mer + 3 variants), under the same eligibility rule — see "Locus eligibility" below, and "Heterozygosity, ploidy, and consensus-assembly inputs" for why the copy-number constraint matters and what it costs.

Canonical invariance: for odd k, the central position maps to itself under reverse-complement (m -> k-1-m = m, base complemented). A transition maps to a transition, a transversion to a transversion — the transition/transversion split is well-defined in canonical space.

Definitions: family, and the canonical form of a family

Family. The family of a k-mer x is the set of (up to) 4 k-mers sharing x's 2m flanking bases, differing only at the central base m. Membership is a property of the flank pattern, not of x itself: any of the 4 possible central substitutions belongs to the same family.

central_canonical_neighbors() (obikseq, CanonicalKmerOf::central_canonical_neighbors) generates all 4 members from any one of them (observed or not), each independently canonicalised (.canonical(), i.e. min(kmer, revcomp(kmer))). This independent canonicalisation is necessary because a central substitution can flip which orientation is lexicographically smaller — two members of the same family can end up canonicalised in different orientations. Despite that, the set of 4 resulting canonical k-mers is invariant: calling central_canonical_neighbors() on any member of a family — present in the index or not — yields the same 4 values. This is relied upon throughout the rest of this document.

Canonical form of a family. Because orientation can differ member to member, "which of the 4 is the reference" cannot be defined relative to whichever member happened to be visited first, nor relative to the minorant (see below) — both are data-dependent (they depend on what is actually observed), so using either as the reference would make the reference itself vary depending on what happens to be present in a given index. Instead: the canonical form of a family is, by definition, the member whose own central base — read in its own already-canonical orientation — is A. This is well-defined for every family, computed purely from the flank pattern, whether or not that specific member (or any member at all) is actually observed anywhere in the index. Concretely: call central_canonical_neighbors() on any member (observed or not) to get the family's 4 canonical forms; the one among them whose own centre nucleotide is A is the family's canonical form. The other 3 (C, G, T) are labelled relative to that fixed reference, not relative to the calling member's own orientation.

Consequence for the minorant. With this fixed A-referenced labelling, minorant (the smallest raw encoding among the family's observed members, introduced further below) becomes directly computable rather than needing to be tracked as extra state: regenerate the family's 4 canonical forms from any member's own k-mer (cheap, no lookup), compare the raw encodings of whichever are marked present, and take the smallest. No separate stored bit is required — see Step 2b below, where this replaces the earlier minorant-bit design.

Locus eligibility: raw definition vs. paralogy filter

For each k-mer x observed in genome A (source, one MPHF slot; the 3 central-position variants generated as in the sweep below): check whether A's locus (flanks fixed) is resolvable in genome B under one of the 4 central forms.

Raw / no model. The locus counts in the denominator iff at least one of the 4 forms is present in B; it counts in the numerator iff the form found in B differs from A's own. No constraint on A's or B's own copy number at this locus. Open question, not resolved: what if more than one of the 4 forms is present in B simultaneously (ambiguous target — count once arbitrarily, count all, or drop)? The stringent filter below sidesteps the question by construction rather than answering it.

Stringent / paralogy-aware. The locus counts only if exactly one of the 4 forms is present in A and exactly one is present in B (count == 1 at that slot too, when a count index is available, to also exclude same-allele duplicates that presence alone cannot see). This drops the raw definition's ambiguous-B case automatically, at the cost of also dropping heterozygous sites indiscriminately alongside true duplications (see "Heterozygosity, ploidy, and consensus-assembly inputs" below).

Rejected: parsimony-based multiset pairing for multiplicity > 1. Rather than dropping ambiguous loci, pair identical alleles between A and B first (0-mutation explanation preferred), then take min(unmatched_A, unmatched_B) as inferred SNP pairs. Rejected on two grounds: (1) circularity — selecting pairs by minimal apparent divergence, then measuring divergence on those same pairs, deflates the estimate by construction, not a neutral heuristic; (2) the discriminating signal is a single base among 4 possible values, and the flanks are already guaranteed identical for every candidate by construction (that is how the locus was selected) — no information remains in a k-mer window to tell which copy in B truly corresponds to which copy in A once multiplicity > 1 on either side. Any pairing rule invents a correspondence the data cannot support. Multiplicity > 1 is treated as non-identifiable, not as a puzzle to solve with a heuristic.

Multi-genome framing: family as pseudo-alignment column

Idea. Instead of resolving locus eligibility and correspondence one genome pair at a time, treat a family as a column of a pseudo multiple alignment across all genomes simultaneously: for each family, each genome has either a net single-copy state (A/C/G/T, when the genome carries exactly one of the 4 forms) or "missing" (?, multi-copy or absent). Flank conservation (the 2m bases fixed by construction) supplies positional homology for free — the same role a real MSA would play, without alignment software, gap penalties, or progressive-alignment approximations. Stacking one such column per family, genomes as rows, produces a genuine SNP pseudo-alignment matrix, not just a bag of pairwise distances.

Precedent. This is the same principle behind reference-free k-mer-based phylogenomics tools — SKA (Split K-mer Analysis, Harris 2018) and kSNP: split the k-mer around a variable center, use flank identity to call homologous columns across arbitrarily many genomes with no reference and no MSA step, then feed the resulting pseudo-alignment to standard phylogenetic tools. Landing on the same design independently is a good sign, not a coincidence.

Resolves the pairwise-correspondence problem, properly. The "Rejected: parsimony-based multiset pairing" case above failed because, with only two genomes' cardinalities to look at, there is no external constraint to justify picking one correspondence between leftover alleles over another — min(a,b) is a lower bound dressed up as a point estimate (see the follow-up discussion on Felsenstein-style parsimony inconsistency: minimum-event explanations are systematically biased low whenever homoplasy/multiplicity is real, not noise-cancelling). With N genomes and many families jointly, the same question can be answered the way real phylogenetics answers it: ancestral state reconstruction / ML mapping over a tree estimated from the whole column set. The tree supplies the missing constraint that two isolated columns cannot — this is the principled way out, not a heuristic replacement for one.

Relation to what's already implemented. KmerIndex::raw_snp_distance already computes, internally, per family, exactly this row — single_form: Vec<Option<u8>>, one entry per genome, None where ambiguous/absent — before immediately collapsing it into pairwise snp[i,j]/shared[i,j] tallies. The pivot this section proposes is small at the implementation level: stop collapsing early, and surface the per-family row as a first-class artifact (a families x genomes matrix). Pairwise raw p-distance becomes one projection of that matrix (what's computed today), not the primary object; downstream, the matrix itself could feed real phylogenetic tools (parsimony/ ML, e.g. RAxML/IQ-TREE-style) instead of only NJ/UPGMA on a homemade pairwise-distance matrix.

Caveat: column completeness shrinks with N. The probability that a family's flanks stay intact simultaneously across all N genomes decays with N (same ascertainment-bias mechanism as Bias 1 above, compounded over more genomes) — fully-resolved columns (no ? anywhere) become rare as more genomes are added. Same missing-data situation any real multi-species alignment faces, and phylogenetic tools already handle it well; the practical 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 phylo --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.

A concrete Sankoff cost matrix for the 16-state alphabet

Why c_gl and c_ctx are the same constant, not two. An earlier version of this design used two independent parameters: c_gl for an ordinary gain/loss step between two nonempty states (e.g. {A,C}->{A}), and c_ctx only for a state collapsing all the way to ∅. That distinction doesn't survive contact with what the index actually observes. Presence/absence tracking never sees "the flanks" separately from "the centre" — it sees whole, distinct, homologous 31-mers (two 31-mers sharing >=30 bases are homologous by construction, though orthology vs. paralogy is undecidable from that alone — settled earlier, see "Locus eligibility"). Losing the A form of a family while C remains ({A,C}->{C}) and losing the last remaining form ({A}->∅) are the same kind of event: a specific, complete 31-mer that used to be observed no longer is, because at least one of its 31 positions mutated. There is no separate "gain/loss of a still-recognised allele" mechanism distinct from "loss of context" — both are exactly the event c_ctx(p) (below) already computes the expected cost of. So: one constant, c_ctx, used everywhere a member is gained or lost — X -> ∅ included, and (see next point) not even as a special case there.

Set-edit-distance formula. For two states X, Y ⊆ {A,C,G,T}, split into seulement_X = X \ Y (size a) and seulement_Y = Y \ X (size b). Elements present in both cost nothing. Pair up to min(a,b) of the remaining elements as substitutions (cheaper than treating them as an unrelated loss + gain whenever c_sub < 2*c_ctx, which any sane parameter choice satisfies); whatever is left over after pairing is a pure gain/loss:

cost(X,Y) = min(a,b)*c_sub + |a-b|*c_ctx

Worked examples (c_sub = c_ctx = 1): {A}->{C} = 1 (one substitution); {A}->{A,C} = 1 (one gain, no substitution pair available since nothing is only-in-Y that matches an only-in-X element after the shared A is excluded); {A,C}->{G,T} = 2 (two substitution pairs, A/C vs G/T, both same-size sets share nothing); {A,C,G}->{A,C,T} = 1 (G/T is the only mismatched pair, A,C shared).

∅ is not a flat-cost special case — a corrected position, reversing an earlier draft of this design. An earlier version charged cost(X, ∅) = c_ctx flat, independent of |X|, on the grounds that every member of a family shares identical flanking sequence, so one mutation breaking that context should unrecognise all of them at once rather than |X| separate times. That argument doesn't survive comparison with how the rest of the matrix already works: {A,C,T} -> {A} (losing two members, one remaining) is charged 2*c_ctx via the general formula above, with no equivalent "maybe it was one shared event" discount — nothing distinguishes that case from {A,C} -> ∅ (losing two members, none remaining) other than which state happens to be the target. Singling out ∅ for special treatment was arbitrary, not principled: the uncertainty about "one event or several" is identical in both cases, and the model already resolves it uniformly elsewhere by simply counting elements. So ∅ is now an ordinary node in the state graph below like any other, with no patch applied afterwards:

cost(X, ∅) = cost(∅, X) = |X| * c_ctx   (via |X| single-element graph edges)
cost(∅, ∅) = 0

Computing c_ctx. Not guessed — built from the value already derived for the >1 bucket in "Context, detectability, and a 3-way ordinal distance per pair" above:

c_ctx(p) = mean_sub_cost * [(2m*p) / (1 - (1-p)^(2m)) + p]

p is p_hat, the calibrated per-site mutation probability (see calibrate_p_hat / the "Experiment" section below); m = (k-1)/2, the flank length on each side of the centre (k=31 → m=15). Read the bracketed term, well, term by term: (1-p)^(2m) is the probability that none of the 2m flanking positions mutated, so 1-(1-p)^(2m) is the probability that at least one did — i.e. the probability that this specific 31-mer stops being observable at all, which is exactly the event c_ctx prices. 2m*p is the unconditional expected mutation count over those 2m positions. Their ratio is the conditional expectation E[mutations | at least one occurred] — provable directly: for any nonnegative integer random variable X, X · 1{X>=1} = X always (both sides are 0 when X=0, both are X otherwise), so E[X | X>=1] = E[X·1{X>=1}]/P(X>=1) = E[X]/P(X>=1), which is (2m*p)/(1-(1-p)^(2m)) here. The trailing +p adds the (much smaller, first-order, not itself a conditional expectation) marginal contribution of the centre position's own mutation probability.

The mean_sub_cost multiplier — a correction, not part of the original derivation. The bracketed term is a count of expected mutations, not a cost — an earlier version used it directly as c_ctx, implicitly pricing every one of those mutations at a flat 1 regardless of type. That stopped being defensible once substitution costs were calibrated per base-pair category (substitution_costs_from_tally): transitions are markedly more frequent than transversions in this project's own data (roughly 5-6x, e.g. 130k vs. 20-28k observed instances per category in the eubacteria run below) and correspondingly cheaper, so "one mutation" isn't worth a flat unit. mean_sub_cost (mean_substitution_cost) is the empirical average substitution cost, weighted by each category's observed frequency — it converts the bracketed term's expected mutation count into an actual expected cost. Concretely, in the eubacteria run: the bracketed term alone is ~1.17 (close to 1, i.e. context loss is usually attributable to a single mutation, since p_hat is small); mean_sub_cost ~= 1.48 (pulled up from the cheapest transition cost of 1.0 by the substantial minority of transversions in the mix); final c_ctx ~= 1.73 — sensibly between the cheapest transition (~1.0) and the transversion costs (~2.5-2.9), rather than coinciding almost exactly with the cheapest transition purely by construction accident, which is what the unmultiplied version did and which is what exposed this gap in the first place.

Better construction method: shortest path in a small state graph, not the closed-form formula directly. Build a graph on all 16 states, ∅ included, with two edge types — substitution edges between same-cardinality sets differing by one element (weight c_sub, or c_ts/c_tv if split further below), and gain/loss edges between sets whose cardinality differs by one (weight c_ctx, per the unification above; ∅ connects to each singleton state this way, being a strict subset of it) — then define cost(X,Y) as shortest-path distance in that graph, uniformly for every X,Y including ∅, precomputed once (16 nodes, trivial) into a dense 16x16 matrix before feeding it to Sankoff/TNT. Verified equivalent to the closed-form formula above in the uniform-cost case (checked by hand on {A}->{C,G,T}: both give c_sub + 2*c_ctx). The graph construction is not just a reformulation for its own sake: it is the version that generalises correctly once substitution costs stop being uniform (next point) — the closed-form's min(a,b) counting silently assumes any pairing costs the same, which breaks the moment c_sub depends on which two bases are involved.

Substitution refinement — implemented as the full 6-category symmetric matrix, not just Ts/Tv. c_sub was originally going to split into just c_ts (A<->G or C<->T) and c_tv (the other four pairs) — already well-defined in canonical space (see "Canonical invariance" above: a transition maps to a transition, a transversion to a transversion, regardless of orientation). A symmetric cost matrix allows finer resolution than that 2-category split, though: each of the 6 distinct undirected base pairs (AC, AG, AT, CG, CT, GT) can be costed independently, with no grouping at all — SankoffWeights::sub_cost is a full [[f64;4];4] table (SankoffWeights::ts_tv remains as a convenience constructor for the coarser 2-category case, if ever wanted). This turns "pick min(a,b) substitution pairs" into a genuine (tiny, <=4 elements per side, trivially enumerable) minimum-cost bipartite matching problem instead of a plain count — the state-graph shortest-path construction handles this automatically, no separate logic needed. Calibration is not a guess either: KmerIndex::base_pair_tally collects the pooled (centre_i, centre_j) distribution over resolved (SNP) sites — restricted to the same ratio-ceiling-included genome pairs as p_hat's own calibration, for the same saturation-exclusion reason — and substitution_costs_from_tally turns that into cost(a,b) = -ln(rate(a,b)), normalised so the most frequent category costs 1.0 (the standard generalised-parsimony step-weighting heuristic, generalised from 2 categories to 6). A second full pass over the annex is required for this (see base_pair_tally's own docs for why it can't share raw_snp_distance's single pass) — both now share one traversal helper, scan_family_pairs, rather than duplicating the per-family reconciliation logic.

Caveat carried over from the D_F = min(a,b) rejection earlier: this cost matrix is only valid as Sankoff step-cost input, re-evaluated for every branch of every candidate topology during tree search. Reusing cost(leaf_A, leaf_B) directly as a standalone pairwise distance (bypassing the tree) would reintroduce the exact circularity already rejected — the min(a,b) pairing here is a locally-defined edit distance between two states, not a claim about the true evolutionary history between two specific genomes.

Feasibility confirmed: TNT's costs command accepts custom step matrices for multistate characters, so this whole construction (16x16 matrix derived from the state graph, c_ts/c_tv/c_ctx as the three tunable parameters) is directly usable there — no new tooling required before testing it.

Experiment: TNT run on real data (2026-08-11)

Status: validated at genus/family/order scale within Bacteria; not informative across domains with this character type. Exploratory only — run entirely outside the repo (/tmp/tnt_run, TNT installed locally under TNT/), no new Rust code. Kept here as the record of what was learned.

Pipeline. obikmer phylo --snp emits one IUPAC-coded pseudo-alignment row per genome (snp_pseudo_alignment, one column per family with family_size() >= 2). A small external Python script decodes IUPAC back to the 16-state bitmask, applies the set-edit-distance formula above, and emits a complete TNT script (xread matrix + smatrix step-matrix + hold/ mult search). No new Rust code was needed for this pass.

Calibration ("quick option"). Rather than modifying write_raw_snp_distance_csv to emit raw counts, c_ctx was approximated as the unweighted mean of the pairwise snp/(snp+shared) ratios already present in rawsnp.csv, restricted to the relevant taxon subset. Used values: mu = gamma = 10, c_ctx = 18 (this run predates the c_gl/c_ctx unification above — gamma here is what c_gl was called before it turned out to just be c_ctx; ratio c_ctx/mu = 1.8, consistent with the observed pairwise ratios for genomes within Enterobacteriaceae/ Eubacteria, which cluster around 1.5-2 and barely move as the taxon set widens — the context signal is stable, not sensitive to which subset is chosen). The more principled route (raw counts, weighted p_hat, Ts/Tv split) remains a follow-up, not yet done.

Run 1 — Enterobacteriaceae (11 taxa: 4 E. coli, 4 Salmonella enterica, 3 Klebsiella pneumoniae). All three genera recovered as monophyletic. Initial read of the exported (unrooted, TNT/Nexus [&U]) tree as showing a genus-arrangement disagreement with known systematics (Escherichieae: Escherichia+Salmonella sister vs. more distant Klebsielleae) was wrong — diagnosed via force = (taxa); monophyly-constraint test (identical score constrained vs. free ⟹ no real disagreement). Root cause of the misreading: with exactly 3 clades and no outgroup, an unrooted tree has only one possible topology (a single trifurcation) — there is no internal arrangement to get right or wrong. This run cannot test inter-genus relationships at all; it can only test intra-genus monophyly (which held).

Run 2 — Eubacteria (18 taxa: run 1 + Acidobacterium capsulatum, Opitutus terrae, Bacillus subtilis, Shouchella clausii, Wolbachia endosymbiont, Proteus mirabilis, Yersinia ruckeri). The Enterobacteriaceae substructure from run 1 is reproduced identically, now correctly rooted by real outgroups, resolving the tribal arrangement left undetermined in run 1: Escherichia+Salmonella sister, Klebsiella more distant — matching known systematics. Outgroup placement: (Bacillus, Shouchella) sister pair (Firmicutes/Bacillales) splits from all Proteobacteria — a correct phylum-level split; Wolbachia (Alphaproteobacteria) splits from the Gammaproteobacteria block (Proteus, Yersinia, Enterobacteriaceae) — a correct class-level split. Lower confidence: the fine nested order (Proteus, (Yersinia, Enterobacteriaceae)) and the relative position of Acidobacterium vs. Opitutus — plausible, not independently verified against current Enterobacterales family-level literature.

Run 3 — full domain set (20 taxa: run 2 + Candidozyma auris [yeast] and Saccharolobus islandicus [archaeon]). The bacterial clade from run 2 is reproduced unchanged and intact — a real robustness signal, the method does not fragment the ingroup when unrelated deep taxa are added. But the result carries no information on Bacteria/Archaea/Eukarya relationships: with exactly one eukaryote, one archaeon and one bacterial clade, the unrooted tree is again forced into the single 3-clade trifurcation from run 1's caveat — there is no second representative of either outgroup domain to resolve internal arrangement, so nothing about their relative position can be read from the topology (a "ladder" ordering in the exported tree is serialization, not signal). Independently, rawsnp.csv shows why this character type cannot reach further: pairwise ratios involving Candidozyma/Saccharolobus are almost all NA (no central-position family shared at all) or saturated at 1.0 (every shared family differs) — central-position families require literal 31 bp context conservation, which simply does not survive domain-level divergence. Cross-domain placement would need conserved-marker characters (rRNA, ribosomal proteins), not this estimator.

Native --sankoff --tnt/--phyg export (2026-08-12), superseding the external scripts above

The ad hoc Python glue from the previous section is superseded: obikmer distance --sankoff now calibrates the matrix natively (p_hat, 6-category substitution costs, c_ctx weighted by mean_sub_cost — see the worked example above) and --tnt/--phyg each write a ready-to-run script from it, no external script needed. ∅ is an ordinary 16th state throughout (never -), specifically to avoid gap-semantics confusion in downstream tools — see "A concrete Sankoff cost matrix" above for why.

TNT (--tnt). write_sankoff_tnt (obikmer/src/cmd/phylo/mod.rs) recodes to TNT's own 0-9A-F xread alphabet (its default reader rejects the wider IUPAC set otherwise), scales and rounds costs to integers (smatrix/cost reject decimals), then re-runs integer Floyd-Warshall on the rounded matrix (scaled_metric_matrix) — independently rounding each cell of an already-metric real-valued matrix can break the triangle inequality (e.g. two real costs of 1.734 round to 173 each, summing to 346, while their real sum 3.468 rounds to 347), which TNT otherwise silently "fixes" itself with an unreproducible correction. Verified against the real 20-genome benchmark index: zero triangle-inequality violations after the fix, TNT loads the file without its "triangle inequality violated... Fixed" warning.

Two syntax facts worth recording because they're wrong in intuitive guesses and contradicted actual TNT behavior when tested: TNT's plain command stream has no comment syntax of its own — /* */ and [ ] only work inside the (separately-enabled, off by default) macro scripting language, and error with "No command!" otherwise. The working substitute is quote TEXT ; (prints the text, doesn't affect parsing) — but the text itself can't contain a literal ; (TNT's universal terminator); the manual's own escape (.,) exists but the script here just avoids semicolons in the text instead.

The default search command embedded in the script is mult (traditional: random addition sequences + TBR), not xmult (New Technology search: ratchet/drift/tree-fusion). xmult with TNT's default mxram (16 MB, must be set before xread if changed) ran out of RAM on the real 908k-character dataset ("xmult - out of ram"); mult does not, matching what had already been validated by hand outside this session.

PhyG (--phyg). write_sankoff_phyg writes a tcm: custom-alphabet matrix (same scale+round+metric-closure treatment as TNT) and reuses --sankoff's own _sankoff.fasta as-is via prefasta: — PhyG's tcm: alphabet is read from the matrix file's own first line, so (unlike TNT) no recoding is needed. PhyG auto-adds its own indel/gap state as an (n+1)-th row/column of the tcm; inert here since the alignment encodes absence as 0, never -.

report("file", newick, overwrite) — exactly as shown in PhyG's own manual — triggers Unrecognized/missing report option ... defaulting to 'graphs' on the locally installed binary (1.3, commit 3c1a1fa); the working form adds graphs explicitly: report("file", graphs, newick, overwrite). Manual/binary mismatches like this (also true of criterion: — the binary accepts parsimony/ml/pmdl, the manual instead documents mapa/ncm/parsimony/pmdl/si) mean command syntax against this PhyG build should be verified empirically, not trusted from the PDF alone. instances:N (not a separate CPU flag) is what actually parallelises the search across cores — PhyG uses all physical cores by default but only across as many instances as are running, so raise it to the physical core count to use them all (the CLI-level +RTS -NX -RTS flag also exists but controls something else: capping/limiting cores, not requesting more).

Both scripts share one --sankoff-cost-scale (default 100), not two separate flags — they scale the same calibrated matrix for the same reason (integer-only cost commands) and no PhyG-specific accumulator-width constraint was ever found to justify a different default from TNT's (TNT: hinted 32-bit accumulators in its own manual; PhyG: no such hint found — Haskell's native Int is typically 64-bit).

Open problem: PhyG reports all branch lengths as 0.0. The graph-level parsimony cost is correct (3.3286×10⁸ on the real dataset, consistent with TNT's 328574911 on the same calibrated matrix), but every individual edge in the exported Newick shows :0.0, with the total cost only ever shown as a whole-tree annotation ([3.32860377e8]). Not fixed, not fully diagnosed — PhyG's manual describes per-edge branch length as computed by ancestral-state (HTU) backtracking, well documented for sequence/standard character types, but nothing found (the term "Sankoff" doesn't even appear in the manual) confirming this backtracking is wired up for a custom tcm: matrix character. Switching criterion: to a likelihood-family option (ml, or the manual's mapa/ncm/si) was considered as a possible fix but is very unlikely to be one: those criteria are information-theoretic reparametrisations of the same step-counting machinery as parsimony (ncm in particular is known in the literature to be numerically equivalent to weighted parsimony), not classical continuous-time-Markov ML with a real rate matrix — so they wouldn't change how branch length is attributed per edge either.

Export format note (not a bug in the generator). Neither script's .tre output opens in PearTree (FigTree's successor) via File > Open — association/Launch-Services quirks were ruled out (the file was opened directly through the app, not by double-click). Likely cause, not yet confirmed: TNT's export is a minimal NEXUS begin trees; block with no preceding Taxa block and bare numeric (untranslated) leaf labels; PhyG's is multiple raw Newick trees concatenated with no NEXUS wrapper at all plus a trailing [cost] bracket tag after the root label. Both differ from a "normal" single-tree, fully-declared NEXUS file; this is PhyG/TNT's own export format, not something write_sankoff_tnt/write_sankoff_phyg could fix without post-processing the other program's output after the fact.

Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)

Decided: parsimony (the whole --sankoff/--tnt/--phyg pipeline above) is a stopgap, not the destination. The goal is maximum likelihood with real, calibrated branch lengths (expected substitutions/site), which parsimony step-counts were never going to give directly (see the open "branch lengths are 0.0" problem above — even if fixed, TNT/PhyG-style parsimony branch length is a step count, not a continuous ML estimate).

Model choices, settled:

• The exchangeability R is symmetric; the rate matrix Q is not. (Superseded an earlier, wrong framing here that treated "symmetric model" as one thing — see the resolution below on R vs Q vs π for the full reasoning.) R(a,b) = R(b,a) because BasePairTally never captured direction — a fact about the data, not a modelling choice. Q(a,b) = R(a,b)·π_b is asymmetric whenever the real state frequencies π are (which they are, empirically) — biology drives this via π, not via R.
• ∅ stays an ordinary 16th state, as already established for TNT/PhyG — same reasoning applies to any ML tool: encode as a real alphabet symbol, never as -/gap, or the RAxML-era failure (empty set silently treated as missing data) repeats.

Stationary frequencies for the 16 states — resolved (2026-08-12). A CTMC needs a rate matrix Q, generally asymmetric. Q(i,j) = R(i,j) · π_j, where R (exchangeability) is symmetric and π (stationary frequencies) need not be — this factoring is what makes Q reversible (satisfies detailed balance, π_i·Q(i,j) = π_j·Q(j,i)) for any π, not just uniform, as long as R is symmetric. Two separate, both-easy quantities, not one hard inverse problem:

• R is already calibrated: sub_cost (-ln(observed rate) per pair, from BasePairTally) is R up to a log transform — recover it as R(a,b) = exp(-sub_cost(a,b)). Symmetric by construction, because the tally itself never distinguished direction (unordered-pair counts only) — not a modelling choice, a fact about what the data can say.
• π is a direct count: empirical marginal frequency of each of the 16 states across the whole alignment (same kind of scan already used to confirm N occurs 1383 times in the real 20-genome benchmark). With ~908k sites × 20 genomes, the counts are large enough that this is precise on its own — no need to spend ML degrees of freedom re-estimating it via IQ-TREE's +FO.

Checked and ruled out along the way: IQ-TREE's +F (empirical, "compute from the alignment") does not work as a shortcut for this — for a custom-file morphology model, readParameters always requires the file's own frequency line unconditionally; omitting it and passing +F instead just fails (ERROR: State frequencies could not be read). π has to be computed by obikmer and written into the file, not left to IQ-TREE.

Net effect: no free-rate ML estimation needed for this piece at all (the mistaken assumption that motivated most of this discussion — that building an asymmetric Q a priori would require solving a linear system from Q itself — doesn't apply, because R, the only piece that's genuinely hard to get directionally, is symmetric and already in hand).

Candidate tool: IQ-TREE, because it supports user-defined multistate models (unlike RAxML's MULTI data type, which is limited to the equal-rate Mk model and can't take a custom rate matrix at all — a genuine tool limitation, not a gap-symbol encoding problem this time). IQ-TREE 3 (3.0.1) is now installed locally (Homebrew, iqtree3).

IQ-TREE custom-model format — verified empirically against the local binary (2026-08-12). The web docs' -mdef NEXUS begin models; frequency NAME = ...; model NAME = ...; end; mechanism (initially assumed to apply directly, see history below) turned out to be for named components used inside MIX{...}/FMIX{...} mixture models only — it does not apply to a single, non-mixture custom morphology matrix, and using it that way fails (ERROR: File not found <name> — traced in IQ-TREE 3's own source, model/modelmorphology.cpp: any -m string that isn't MK/ORDERED/GTR/GTRX is passed straight to ModelMarkov::readParameters(), which opens it as a literal file path, never consulting the -mdef models block at all for this data type).

The confirmed working recipe (built a tiny 5-taxon/3-state toy dataset and rate file, ran it end to end with iqtree3, got a real ML tree with non-zero branch lengths and an optimized log-likelihood — ground truth, not documentation):

• No -mdef needed. Write one plain file (any name) containing, as whitespace/newline-separated numbers, in order: the lower-triangular rate matrix (N(N-1)/2 values, PAML row-major order — for 16 states, 120 values, the same count and layout already produced for TNT's smatrix), immediately followed by the N state frequencies on the same stream (no header, no blank line required — confirmed by reading ModelMorphology::readRates/ModelMarkov::readStateFreq directly, which just pull tokens off the stream in sequence).
• Invoke with -m <path-to-that-file>+ASC (+ASC for the no-constant-site correction, as before). An explicit +F{f1,...,fN} on the command line overrides the file's own frequency line if given (confirmed in ModelMorphology::init) — useful once real calibrated stationary frequencies exist, a placeholder equal-frequency line works meanwhile (the still-open gap noted above).
• --seqtype MORPH (alphabet 0-9,A-V, ≤31 states) — reuse the same 0-9A-F recoding already built for TNT (TNT_STATE_SYMBOL).

Risk, confirmed, precisely characterised, and resolved by design (2026-08-12). --seqtype MORPH{16} does not force the state count for real ML analysis — tested directly (--seqtype MORPH{4} on the 3-symbol toy alignment gave the byte-for-byte identical 3-state result as no {4} at all) and confirmed in source: the value it sets (params.alisim_num_states_morph, utils/tools.cpp) is consumed only by the --alisim simulator; the main analysis path always calls getDataBlockMorphStates/an equivalent scan (alignment/alignment.cpp), for both FASTA/PHYLIP and NEXUS input (a NEXUS symbols= declaration doesn't change this either — checked, same code path). No CLI flag or NEXUS declaration overrides it.

The precise rule (from getDataBlockMorphStates, alignment.cpp:1058): N = one plus the highest state ordinal actually observed anywhere in the alignment, ordinal being the symbol's position in IQ-TREE's own fixed table "0123456789ABCDEFGHIJKLMNOPQRSTUV" — not a count of distinct symbols seen. So the risk is narrower than "any missing symbol breaks it": concretely, it's whether the symbol mapped to state index 15 (F in the 0-9A-F recoding already used for TNT, i.e. our N = "all four bases ambiguous") occurs at least once anywhere in the real alignment — if it does, N correctly comes out to 16 regardless of which lower-index symbols (including 0/∅) are rare or absent; if it doesn't, N silently undercounts and misaligns every value in a 16-entry rate/frequency file, with no error to catch it. Checked against the real biological alignment (/tmp/msg_test/eub_sankoff.fasta, the 20-genome benchmark index): the symbol N (IUPAC "all four bases ambiguous," state index 15) occurs 1383 times across 13 of the 20 sequences — present, so this specific real dataset is not at risk. Still worth a real, general presence check inside obikmer before this is wired in, rather than assuming every future dataset will have N too (nothing in IQ-TREE would catch it if not).

Resolution: subset + compact-renumber, not rely on all 16 appearing. Since IQ-TREE always infers N from the alignment's own content and nothing overrides that, the fix is to make the file obikmer writes match that inference by construction, for every run, rather than hope the 16th (or any particular) state happens to occur:

1. Scan the real alignment for which of the 16 canonical states actually occur anywhere (not per-column — anywhere in the whole alignment).
2. Renumber the occurring states to a compact, consecutive 0..k-1 range, preserving their relative order (the original bitmask/STATE_ SYMBOL ordering) — not just filtering, since a gap in the ordinal sequence (e.g. keeping states {0,1,2,4} numbered as-is instead of {0,1,2,3}) reproduces the exact same "highest observed ordinal" miscount this was meant to fix.
3. Recode the alignment itself with this new compact k-symbol alphabet (same recoding mechanism already used for TNT's 0-9A-F, just over a possibly-smaller symbol set).
4. Extract the matching k×k submatrix (rows/columns for the kept states only) from the full calibrated 16×16 cost matrix, in the same lower-triangular order the rate-matrix file needs — and, later, the matching k-length subset of stationary frequencies once those are calibrated (still the open gap noted earlier in this section).

Consequence, and why nothing is lost: a state that never occurs in a given alignment can, by definition, never contribute a transition to score in that same alignment — dropping it from that run's matrix costs nothing. The subset (and therefore k) can differ from one dataset/run to the next; this has to be done freshly per alignment, not computed once and reused.

(Superseded reasoning, kept for the record: the -mdef NEXUS route below was the original plan, based on IQ-TREE's own web documentation for protein mixture models, before the local install allowed testing it — GTRX combined with a -mdef-referenced custom model, +Fname frequency reference. Both pieces exist and parse without error individually, but GTRX/GTR are IQ-TREE's own fixed built-in equal-structure multistate model, not a hook for an arbitrary custom matrix; a custom matrix is a file path in -m directly, no -mdef or GTRX involved.)

Source for the empirical findings above: model/modelmorphology.cpp and model/modelmarkov.cpp in the local iqtree/iqtree3 source (cloned to inspect the exact parsing logic after documentation didn't resolve the +Fname reference error) — more reliable here than the PDF/web manual, which (like TNT/PhyG) doesn't always match this specific binary. Original (partially superseded) sources: Substitution Models, Complex Models.

Relation to the existing calibration. sub_cost[a][b] = -ln(rate) (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. 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/phylo/iqtree.rs) implements exactly the R = exp(-cost) / empirical-π design above: writes <prefix>_iqtree.model (lower-triangular R, PAML order, then π) and <prefix>_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, <prefix>_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 phylo) 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 -lns) — 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.

Implemented (2026-08-12): pairwise_cost_matrix replaces build_cost_matrix entirely

New module obikindex/src/cardcomp.rs, replacing sankoff::build_cost_matrix and the c_ctx/SankoffWeights/PHatEstimate/calibrate_p_hat/ c_ctx_from_p_hat/mean_substitution_cost/substitution_costs_from_tally machinery it depended on outright — not kept in parallel as a fallback (all now unreferenced outside their own tests; sankoff.rs itself is a pending removal, not yet done).

New primitives, obikindex/src/siblings.rs:

• BasePairTally gained a same: [u64; 4] field (diagonal — "both genomes at the same single base", pooled from base_pair_tally's existing traversal, extended to also tally the bi == bj case it previously discarded).
• CardinalityTally { counts: [[u64; 5]; 5] } and KmerIndex::cardinality_tally, a new traversal (same shape as snp_pseudo_alignment's — needs full per-genome presence masks, not scan_family_pairs's single-resolved-form view, since cardinality 2-4 is exactly the signal being tallied, not noise to drop). Same saturated/no-data pair exclusion as base_pair_tally. Restricted to variable families (family_size() >= 2), matching snp_pseudo_alignment's own scope. Verified against the 20-genome benchmark: counts match the earlier hand-rolled Python analysis exactly (e.g. c=0/c=0: 117,158,166, c=0/c=1: 13,707,223 — the same numbers this whole investigation started from).

cardcomp.rs:

• cardinality_transition_probs/composition_transition_probs: row- normalise the two tallies into proper transition probability matrices, diagonal included ("stay the same" is a real, calibrated outcome).
• pairwise_cost_matrix: for every pair of the 16 states, shared = A∩B contributes ∏ P_composition(x→x); lost = A\B, gained = B\A are parsimony-paired (best_pairing_cost, brute-force over the ≤4! injections — small enough that hand-rolling beats a dependency) into substitution events on P_composition, minimising total -ln; the cardinality-difference leftover is priced once via P_cardinality(|A|→|B|), never chained. Row-normalised, -ln'd, then symmetrised: cost_sym(A,B) = (cost(A,B)+cost(B,A))/2 — equivalent to taking the geometric mean of the two raw probabilities (-ln(√(P(A,B)·P(B,A))) = (-ln P(A,B) - ln P(B,A))/2), not their arithmetic mean. Required, not just convenient for IQ-TREE's lower-triangular file format: Sankoff parsimony's score is independent of where an unrooted tree (what TNT/PhyG actually search over) gets rooted only if the cost matrix is symmetric — the discrete-parsimony analogue of CTMC reversibility, established by direct reasoning, not assumed. Bonus of the same decision: 120 free parameters instead of the 240 a fully asymmetric matrix would need.

Verified on the 20-genome benchmark: resulting matrix symmetric (checked numerically, zero asymmetric cells), zero diagonal, no NaN/Inf. --tnt output still loads into TNT with no triangle-inequality warning (scaled_metric_matrix's rounding-metric-closure step still needed and still applied — nothing in the new construction guarantees the rounded integer matrix stays a metric, even though the real-valued one is exact by construction here, unlike the old Floyd-Warshall-closed matrix which needed it for a different reason). IQ-TREE loads the new model file and reports the same π as before (only R changed).

Two consistency bugs found and fixed post-implementation (2026-08-13)

--exclude-genome didn't drop columns that become monomorphic once the excluded genome(s) are gone. snp_pseudo_alignment's "variable family" test (family_size() >= 2) is a property of the annex computed over every genome in the index — unaffected by the CLI-level exclusion, which only dropped the excluded genome's row. A family variable only because of the excluded genome stayed in the alignment as a now-constant column — silently wrong data for TNT/PhyG, and a hard failure for IQ-TREE's +ASC (verified: excluding 2 taxa on the benchmark left 116,351 such columns — matches the manual +ASC failures hit earlier in this same investigation, before --exclude-genome existed). Fixed in drop_excluded (obikmer/src/cmd/phylo/mod.rs): after dropping excluded rows, re-scan each column among the surviving sequences and drop any that are now constant. Verified: 908,723 → 792,372 sites after excluding 2 taxa, zero monomorphic columns remain, π recomputed from the corrected alignment matches an independent recount exactly. The compact-alphabet renumbering (iqtree::compact_alphabet) needed no equivalent fix — it already recomputes which of the 16 states occur fresh on every call, from whatever alignment it's actually handed, so a symbol disappearing (e.g. excluding every genome that carries N) is already handled correctly; verified directly (excluded 13 genomes to force N out: "15 of 16 states" reported, correctly-shaped model file).

cardinality_tally's family_size() >= 2 filter looked inconsistent with base_pair_tally — removing it was tried, and was wrong; reverted. cardinality_tally (modelled after snp_pseudo_alignment) had the filter; base_pair_tally didn't (it visits every family via scan_family_pairs unconditionally, folding fully-invariant loci into its own same diagonal). Read as cardinality_tally under-counting its diagonal relative to base_pair_tally, and — independently — as another angle on the --exclude-genome drift (family_size() being global-only meant a family kept here post-exclusion could differ from what the now-correctly-filtered alignment kept). First fix tried: drop cardinality_tally's filter entirely, matching base_pair_tally's whole-annex scope.

That fix was empirically wrong, confirmed by a real IQ-TREE run, not just a hunch. Log-likelihood dropped from the earlier correct run's -8,364,671/-8,371,082 to -9,170,228 (worse fit, not better), with repeated NNI search needs unusual large number of steps (20) to converge! warnings — and the completed run's total tree length came out at 67.644, roughly 30× the earlier correct runs' ~2.0, i.e. branches blowing up/saturating. Root cause, only clear in hindsight: +ASC ("ascertainment bias correction") exists specifically because the likelihood only ever sees variable sites — the alignment fed to IQ-TREE, by construction, contains not one invariant column. Calibrating R from a population overwhelmingly dominated by genome-wide invariant background (family_size()<2 loci outnumber the ~908k variable ones by orders of magnitude) describes a completely different population than the one +ASC and the alignment actually model — the "consistency" argument for matching base_pair_tally's scope was real, but pointed the wrong way: base_pair_tally's own unrestricted same diagonal turned out to have the identical latent bug (only unmasked once its diagonal existed at all, which happened earlier the same day when same was added), not a correct baseline to match cardinality_tally to.

Final fix: restored cardinality_tally's family_size() >= 2 filter, and gave base_pair_tally's same diagonal the equivalent restriction — scan_family_pairs (shared with raw_snp_distance, which legitimately does want fully-invariant families counted as shared) now passes an extra variable: bool (the family's own family_size() >= 2) to its on_pair callback; base_pair_tally only increments same when variable is true, raw_snp_distance's callback ignores the new argument. Both tallies now describe the same variable-families-only population the +ASC-corrected alignment does. Verified: calibration counts back to their original values exactly (c=0/c=0: 117,158,166, matching the pre-regression run bit for bit), and a full IQ-TREE rerun converged normally — log-likelihood -8,389,106.273 (same order as the two earlier correct runs), total tree length 2.044 (was 67.644), no NNI convergence warnings.

Heterozygosity, ploidy, and consensus-assembly inputs

A within-genome multiplicity signal (more than one of the 4 central forms present at a locus) is produced identically by two distinct causes: paralogous duplication and diploid/polyploid heterozygosity. K-mer data alone cannot distinguish them. The one real discriminator is sequencing depth (heterozygous site: total depth of the present forms ~= the genome's single-copy average; duplication: ~2x or more) — but that signal only exists if genome "counts" are raw-read depth (FASTQ input), not occurrence counts in an assembled FASTA, where per-locus depth is not preserved.

Magnitude is taxon- and mating-system-dependent, not universal. Heterozygosity density: mammals ~1 site / 1-1.5 kb (~0.1%); highly outcrossing plants (maize, poplar) reported an order of magnitude higher (~1%); self-fertilising plants (Arabidopsis thaliana) near zero — but with a documented failure mode where segmental duplication masquerades as "pseudo-heterozygosity"; fungi split between haploid vegetative stages (non-issue) and dikaryotic Basidiomycetes, where two long-diverged haploid nuclei coexist without fusing. The estimator's target use case (closely related genomes, k=31) is exactly where the stringent filter above costs the least for low-heterozygosity taxa and the most for outcrossing/dikaryotic ones — no universal threshold; this is a scope caveat to document, not a problem to solve generically.

Why assembled-consensus inputs don't make measured distances wrong. Phylogenetic inputs are near-universally assemblies, not raw reads, and assemblers collapse heterozygous sites to one consensus allele per position — effectively an arbitrary, largely uncorrelated-between-assemblies choice at each het site. This does not inject unbounded noise: standard population genetics gives d_xy = d_a + (pi_A + pi_B)/2 — the expected pairwise difference between a random allele of population A and a random allele of population B equals the net (fixed) divergence d_a plus the average of the two populations' own within-population diversity pi. Consensus flattening realises exactly this random-allele draw, so the measured genome-to-genome distance is a d_xy-like quantity, not d_a — inflated by heterozygosity by a well-characterised additive term, not distorted unpredictably. The term is negligible when pi << d_xy (the common case for cross-species comparisons), and becomes material precisely in the two cases already flagged above: very closely related genomes (this estimator's explicit target) and highly heterozygous outcrossing organisms, where pi and d_xy are the same order of magnitude.

Caveat: this assumes the flattening is uncorrelated with the phylogenetic signal — plausible for de novo assembly, not guaranteed for reference-guided assembly biased toward one allele (e.g. the reference's) at each het site, which would turn the noise term into a systematic bias toward the reference lineage. Not evaluated here.

Forward-looking implication, not part of the current design. The multiplicity > 1 signal discarded by the stringent filter is a crude per-genome proxy for pi (under low background paralogy). If a pi_hat per genome were tallied alongside SnpTally, a d_a correction (p_hat - mean(pi_hat_i, pi_hat_j)/2, roughly) could recover an estimate closer to net divergence instead of d_xy — a possible extension, not scoped here.

Sufficient statistic: 4x4 base-pair tally

Tabulating the joint distribution of (center_i, center_j) over conserved-flank loci, per genome pair, is sufficient for every downstream correction:

Estimator	Input	Formula
Raw p-distance	total off-diagonal / total	p = SNP / (SNP + shared)
Jukes-Cantor	p	d = -3/4 * ln(1 - 4p/3)
Kimura 2-parameter	transition rate P, transversion rate Q	d = 1/2 ln(1/(1-2P-Q)) + 1/4 ln(1/(1-2Q))
LogDet/paralinear	full 4x4 + base-composition margins	d ~= -1/4 ln det(F), robust to non-stationary base composition

JC/K2P need only the total and the transition/transversion split (the diagonal collapses to a single "shared" total). LogDet needs the full 4x4, already populated at no extra cost (see Step 1/2 below).

Memory for the 4x4 tally: n^2 * 16 counters. Trivial for the project's genome-scale use case (tens to hundreds of genomes); ~13 GB at n=10^4 — outside scope but worth flagging if n grows.

Biases (properties of the estimator, not defects)

1. Conserved-flank ascertainment bias. Only SNPs with intact 2m-base flanks are visible; window-intact probability decays as (1-p)^{2m}. For k=31 (2m=30): 0.74 at p=1%, 0.21 at p=5%, 0.04 at p=10%. This estimator targets closely related genomes. Under rate heterogeneity across sites (universal in practice), conserved flanks correlate with slow centers, so p_hat underestimates the genome-wide average rate — it specifically estimates the substitution rate of conserved regions. Two distinct factors are at play here, not one: P(centre of a given window is a SNP) = p exactly, independent of k — a direct restatement of the raw per-site rate via the bijective window<->centre-position correspondence (Statistic section above), not a k-dependent quantity. (1-p)^{2m} is the separate, genuinely k-dependent ascertainment factor (are the flanks also intact). The two multiply: P(usable window showing a central SNP) = p * (1-p)^{2m} — e.g. at p=1/31 (~3.2%), k=31: p * (1-p)^30 ~= 0.0323 * 0.374 ~= 1.2%, i.e. about 1 window in 83, not 1 in 31 (which is only the centre-mutated fraction, before requiring intact flanks).
2. Bias toward isolated SNPs. Two SNPs within k of each other disqualify each other's flanks. Hypervariable regions are invisible by construction.
3. Indels are invisible. A frameshift destroys k-mer matches in a block; this channel captures substitutions only. Indel divergence shows up as lost shared k-mers (lower Jaccard/Mash), not as SNP signal.
4. k-dependent specificity. "A k-mer match implies common ancestry" is quantitative. For a 3 Gbp genome, expected random flank-30 collisions (k=31): (3e9)^2 / 4^30 ~= 8 — negligible. At k=21: (3e9)^2 / 4^20 ~= 2e6 — no longer negligible. k=31 is safe; k<=21 is marginal to unreliable for large genomes. The large k that guarantees homology is the same k that shrinks the detectable-divergence window — an inherent tension.

Implementation: avoid materializing a de Bruijn graph

A central-SNP pair is topologically a simple bubble in the colored de Bruijn graph (source/sink k-mer shared, two length-k branches differing only at the midpoint). Classical bubble-calling (Cortex/discoSNP-style) finds these, but requires the graph — nodes plus adjacency for ~10^9 colored k-mers — resident in memory. Rejected: prohibitive RAM for this project's scale.

A naive per-pair generalisation of variant lookup across n genomes (query each non-shared k-mer's 3 central variants against every counterpart genome's index) costs O(n^2 . N . 3) random lookups, with the same k-mer's 3 variants regenerated and requeried once per counterpart genome — pure redundant work. Rejected as the basis for an n-genome design.

Implementation: sequential per-partition sweep (no scratch, no graph)

KmerIndex::distance() already opens every partition's presence_store/ count_store simultaneously, memory-mapped, into one LayeredStore (distance.rs:73-77). "Querying another partition" is therefore not a new I/O pattern to design — it is the same O(1) MPHF+evidence lookup the query command already performs at scale. This lets the SNP tally be computed with no scratch files and no auxiliary graph, by sweeping partitions once each as a source:

1. For source partition p, enumerate its distinct k-mers (one per MPHF slot; each already carries its full multi-genome presence/count vector — no need to explode per (k-mer, genome) occurrence).
2. For each, generate the 3 central-substitution variants and canonicalise each independently (min(kmer, revcomp), exactly as any normal query) — this avoids the orientation edge case a masked-flank grouping would have (a substitution that flips canonical orientation is handled correctly because each variant is canonicalised on its own, not inferred from a fixed-orientation flank key).
3. Compute each variant's target partition q via its minimizer; batch/sort the partition's outgoing variant queries by q for locality.
4. Look up each variant in q's already-mmap'd MPHF+evidence; on a hit, combine the source's presence vector (base a) with the variant's presence vector (base b): for every i carrying a and j carrying b, tally[i,j][a,b] += 1.

Deduplication needs no persisted state. Sweeping partitions in a fixed order p = 0, 1, ..., P-1 and only acting on a variant when its target partition q >= p guarantees each unordered SNP pair is counted exactly once: a pair with q < p was already resolved earlier, when q was itself the source partition and p (being >= q) was a valid forward target. No cross-partition flag array is needed — the sweep order is the deduplication rule. Within the same partition (q == p), a lightweight transient tie-break suffices: either a #slots(p)-bit scratch flag reset per partition, or simply comparing the two k-mers' raw u64 encodings and only counting when kmer_source < kmer_variant — no storage at all.

This is a distinct computation stage, not a partial_* in the existing additive-by-partition sense: step 3-4 read across partition boundaries by construction, unlike the row-local partial_jaccard/partial_threshold_jaccard primitives. But it requires no new index files, permanent or scratch: unitigs.bin, mphf.bin, evidence.bin, and the presence/count columns are read as-is, and the only extra memory is the current partition's small outgoing-query batch (#kmers(p) * 3, released once p is done) plus the persistent tally accumulator (n^2 * 16 counters, see above).

Outer loop (over source partitions p) must stay sequential. Two independent reasons, not just one: (a) memory — the bounded-footprint claim above only holds with one partition's outgoing-query batch in flight; running T source partitions concurrently multiplies that batch by T, exactly the blowup the design avoids; (b) correctness — the q >= p deduplication rule requires partitions to be claimed as sources in a fixed order; running p1 < p2 concurrently gives no guarantee p1 has finished claiming its q >= p1 targets before p2 starts claiming its own, breaking the "counted exactly once" property.

Inner loop (target-partition lookups for a fixed p) parallelises safely. Each lookup is an O(1) read against an already-mmap'd structure, independent of the others, with no growing allocation — no memory blowup, no ordering dependency between different q. The only shared mutable state is tally; give each worker thread a thread-local partial tally (fixed n^2 * 16 size, independent of partition size) and merge into the global tally once p's inner loop completes — the same reduce-then-merge pattern Rayon already uses elsewhere in this codebase to open partitions in parallel. Extra memory: #threads * n^2 * 16, negligible (512 MB at n1000, 32 threads) and unrelated to partition size.

Cost: 3 * N_distinct MPHF lookups total across the whole index (each partition swept once as source) — the same order of magnitude and the same operation as running query over the index's entire k-mer content against itself, three times. This is the tool's already-optimized regime, not a new I/O profile to validate.

Cheaper: subsampling

Since the target is a ratio, restricting the source-partition sweep to a bottom-s hash sketch (only enumerate k-mers with hash < threshold as sources) divides the lookup count by the sampling factor without biasing p_hat. Mash-like tradeoff: rate estimated from a sample, not the full k-mer set.

Recommendation

Sequential per-partition sweep (Route D): reuse the already-mmap'd per-partition MPHF/evidence/presence structures for O(1) variant lookups, dedup via the fixed sweep-order rule (q >= p, plus an in-partition tie-break), no scratch files, no graph materialisation. Both the SNP (off-diagonal) and shared (diagonal) counts are accumulated by this same sweep, under the locus-eligibility rule chosen (raw or paralogy-filtered) — not reused from the general-purpose shared_kmers matrix, whose raw-identity definition does not apply the same copy-number constraint. Distances (p, JC, K2P, LogDet) as finalisations of the resulting 4x4 tally, mirroring the partial_* -> *_dist_matrix pattern used for Jaccard/Mash/Bray-Curtis/etc.

Detailed implementation plan

Grounded in the current codebase. File/type references are anchors, not prescriptions; adjust to reality when implementing.

Step 0 — new low-level primitives (obikseq)

Two helpers do not yet exist and are prerequisites:

1. Central neighbours. CanonicalKmerOf<L> already exposes left_canonical_neighbors() / right_canonical_neighbors() (obikseq/src/kmer.rs), each returning the 4 canonicalised neighbours at an end position. Add central_canonical_neighbors() returning the 4 variants at position m = (k-1)/2 (each independently canonicalised via .canonical()). The 3 that differ from the source are the query variants; skip the identity. Building on nucleotide(i) / the raw 2-bit layout keeps it O(1).
2. Lone-k-mer minimiser. Routing a synthetic variant to its partition needs its minimiser, but RollingStat (obiskbuilder/src/rolling_stat.rs) only computes minimisers incrementally along a sequence. Add a standalone minimizer(kmer) -> Minimizer that scans the k-m+1 m-mer windows (PackedSeq::mmer, obikseq/src/packed_seq.rs), canonicalises each, and takes the min by seq_hash() — the same selection RollingStat performs, evaluated once. Partition index is then (minimizer.seq_hash() & (n_partitions - 1)) as usize, exactly as QueryBatch::from_records (obikmer/src/cmd/query.rs:142); n_partitions is a power of two so the mask is valid.

Step 1 — the tally accumulator (obikindex)

A SnpTally holding, per genome pair, the 4x4 joint count of central bases: n * n * 4 * 4 u64 (or a packed lower-triangular form since it is symmetric). Provide merge(&mut self, other: &SnpTally) for the thread-local reduce, and accessors yielding, per pair (i,j): total off-diagonal (SNP), diagonal (shared, i.e. p_hat's denominator minus SNP), transition count P, transversion count Q. The diagonal is always populated — it is not an optional LogDet-only extra, since p_hat's denominator is no longer sourced from the external shared_kmers matrix (see "Locus eligibility" and "Statistic and correspondence with shared" above): the source k-mer's own presence/count vector, already in hand when it is enumerated, supplies the diagonal entry directly, at no extra lookup cost.

Step 2 — the sweep (obikindex, new snp.rs)

Mirror distance.rs: open the presence or count store per partition. But instead of a per-partition partial_*, run the sequential source sweep:

for p in 0..n_partitions:                       # OUTER — sequential
    open source partition p's layers (QueryLayer-style, obikpartitionner)
    enumerate distinct canonical k-mers of p (one per MPHF slot) with their
        presence/count vectors                   # column-major, as query stage 2
    par_iter over these source k-mers:           # INNER — rayon, thread-local tally
        apply eligibility rule to the source's own vector (raw: none;         # diagonal
            stringent: exactly one of the 4 forms present in each genome)     # gate
        for i in genomes eligible with source base a:
            for j in genomes eligible with source base a:
                thread_tally[i,j][a,a] += 1        # diagonal — no extra lookup
        for each of the 3 central variants:
            q = partition_of(variant)
            if q < p: continue                    # dedup: forward targets only
            if q == p and variant <= source.raw(): continue   # in-partition tie-break
            slot = layers[q].find_slot(variant)   # MphfLayer::find, mmap'd
            if hit:
                vb = variant presence/count vector
                apply eligibility rule to vb (as above)
                for i in eligible genomes with source base a:
                    for j in eligible genomes with variant base b:
                        thread_tally[i,j][a,b] += 1
    merge thread-local tallies into global SnpTally

The inner lookup is precisely QueryLayer::find_slot + col_value(g, slot) (obikpartitionner/src/query_layer.rs) — reuse or factor out that path rather than reimplementing MPHF access. Enumerating "all distinct k-mers of a partition with their vectors" is the dump/query stage-2 column-major scan already implemented in dump_layer.rs / query_partition_with; factor a reusable iterator if none fits.

presence_threshold applies exactly as elsewhere: a genome "carries base b" iff its count at that slot is >= presence_threshold (trivially >= 1 for presence indexes).

Open problem (unresolved, session end — not yet fully convinced)

The q >= p / tie-break dedup rule in Step 2's pseudocode above is flawed for the stringent (paralogy-filtered) eligibility rule: it only ever brings two family members into view at once (the source and one looked-up variant), never all four simultaneously, and which subset gets compared depends on partition sweep order. "Exactly one of the 4 forms present in genome A" is a whole-family property and cannot be decided correctly from a sequence of pairwise, order-dependent glimpses — the pseudocode above needs revision, not just the eligibility gate bolted onto it as written.

Direction discussed, not yet settled:

1. Every distinct source k-mer looks up all 3 variants unconditionally (drop the q < p skip entirely) so that every observed family member independently gathers all 4 vectors (its own + whichever of the 3 variants exist) at once — a whole-family, order-independent view, computed redundantly once per observed member. Same total lookup order of magnitude as already budgeted (3 * N_distinct), just organised differently (no lookup actually skipped, versus the original rule which skipped roughly half).
2. Tie-break after gathering, not before: only the member whose own canonical encoding is the smallest among the members actually observed (now known, since all were just looked up) writes to SnpTally; the others silently discard their redundant computation. Deterministic, order-independent — as a side effect this also removes the "outer loop must stay sequential" constraint from the cost/parallelism discussion above, since no step depends on partition processing order any more.
3. Proposed optimisation: precompute, once at index build time, a compact global (not per-genome) annex per MPHF slot — the count of other family members observed anywhere in the dataset (0-3). Slots with count 0 (majority under low divergence and few genomes, but see the scaling caveat below) need no cross-lookup at all: eligibility reduces to a local count == 1 check at that single slot, and only slots with count

   = 1 enter the 3-lookup sweep machinery above. Revised (see Step 2b below): minorant status is stored alongside the count after all, on 3 bits rather than 2 — it comes for free from the same lookups needed to count siblings, and storing it lets the sweep discard non-minorant slots without re-fetching anything.

Minorant/sibling-count relationship, worked out precisely. "Minorant" is a one-way implication from sibling count, not an equivalence: 0 siblings => minorant (trivially — with no other observed member, the k-mer is by definition the smallest of the observed set, itself alone), and its contrapositive not minorant => >= 1 sibling. The converse does not hold: being the minorant says nothing about sibling count — a minorant can have 0, 1, 2 or 3 siblings, all with larger encodings than itself. Consequence: this confirms, as a logical necessity rather than a heuristic, that a 0-sibling slot can always write its diagonal contribution with zero ambiguity and no lookup (it is unconditionally its own minorant) — but it gives no shortcut for the >= 1-sibling case, where minorant status still requires the actual comparison of gathered encodings; sibling count alone never determines it.

When to compute the annex, and cache invalidation. Sibling count is a property of the whole set of columns (genomes/groups) currently in the index, not of any single genome — it cannot be computed correctly at mono-genome build time (a family may gain siblings, or its minorant may change, once more genomes are merged in later). Computing it eagerly at every merge would also waste work on intermediate merged states nobody ever queries. Instead: compute it lazily, on first phylo call against a given index, and persist the result alongside that index for subsequent calls — the same lazy-derived-cache pattern PersistentBitMatrix already uses for Columnar -> Packed. This requires no explicit invalidation for merge or filter (obikindex/src/merge.rs, obikmer/src/cmd/filter.rs): both only ever write to a fresh --output directory, never mutate an input index in place, so a re-merged/re-filtered index is simply a new state with no annex yet. select --in-place (select_layer.rs:139-235) is the exception: it aggregates genome columns into groups (Any/All/None/Sum/Min/ Max) by mutating the existing index's files without changing its location. It does not remove k-mer rows, but it can still change eligibility and sibling counts derived from those rows (e.g. a Sum over several single-copy genomes can read as multi-copy at the group level). Because it mutates in place, select --in-place must explicitly invalidate (delete or mark stale) any cached sibling-count annex for that index — the one operation in the current pipeline where this doesn't happen for free.

Not yet convinced this is the right shape, and Step 2's pseudocode above has not been rewritten to match — flagged for the next pass rather than resolved here.

Step 2b — sibling-count / minorant annex (consolidated plan)

Scope: only the precursor annex — not the SNP tally itself, whose Step 2 sweep remains unresolved above. This piece is simpler than the sweep, because it writes to an independent per-slot value, not a shared cross-k-mer accumulator, so it needs no dedup/ownership logic at all at this stage.

Revised annex encoding — 4-bit presence mask, not 3-bit (minorant + count). Superseded after settling the "canonical form of a family" definition above. The 3-bit design (1 minorant bit + 2-bit sibling count, § below, kept for the historical record) had two problems: it discards which variants are present (only how many), so any future consumer (the SNP sweep, or a stats pass — see below) that needs to know which bases exist still has to regenerate and blindly re-query all 3 candidates; and the minorant bit's meaning was tied to whichever member was visited, not to a fixed reference. Storing instead a 4-bit mask — one bit per base (A/C/G/T), set iff that member of the family (labelled relative to the family's fixed canonical form, i.e. the member with A at the centre — see above) is observed anywhere in the index — fixes both:

• Sibling count is derived, not stored: siblings = popcount(mask) - 1.
• Minorant is derived, not stored: regenerate the family's 4 canonical forms from the slot's own k-mer (cheap, no lookup — see above), compare the raw encodings of whichever bits are set in the mask, take the smallest. Erratum (2026-08-14) — this was wrong, kept struck through rather than deleted. "Cheap, no lookup" only accounts for the bit algebra (regenerate 4 forms, compare raw encodings) — true in isolation, but it silently assumed "the slot's own k-mer" is a free fact. It isn't: getting from an MPHF slot index back to the actual k-mer sequence means reconstructing slot_kmer for the whole layer — scan unitigs.bin, mphf.find() every k-mer to place it — an O(distinct k-mers) pass through the MPHF, not a per-slot O(1) lookup. That reconstruction is free only when the caller already needs k-mer identity for something else in the same traversal (e.g. the SNP sweep below, which needs it anyway to generate central_canonical_neighbors()). A caller that wants only the minorant flag pays the full reconstruction for nothing: measured on a real run, a bare family-size histogram (four buckets, otherwise near-instant) spent 71% of wall-clock in MphfLayer::find, all of it solely to answer "is this slot the minorant". Current design: minorant is stored after all — a 5th mask bit (4 presence bits + 1 minorant bit, still fits one byte alongside the presence mask below), written once in the construction pass where the k-mer is already in hand for other reasons (obikindex::siblings::build_layer_sibling_annex), read back for free by every later consumer (FamilyMask::is_minorant). This is genuinely a return to the superseded 3-bit design's core idea (store minorant alongside the count) — the "which variant" blindness that motivated moving away from it is fixed by keeping the full 4-bit presence mask too, not by dropping the stored minorant bit again.
• A future consumer knows exactly which variants to (re-)query — popcount(mask) - 1 lookups instead of always 3, and it knows which 3 (or fewer) to issue, not just how many hits to expect.
• The all-zero value (no base present at all) is still logically unreachable as a real result — the slot's own base is always present in its own family — so it remains available as a free "not yet computed" sentinel, exactly as before.

1. Primitive. Reuse central_canonical_neighbors() from Step 0 unchanged — the 3 canonicalised central-substitution variants of a k-mer (plus the identity, i.e. all 4 members of the family — see "Definitions" above).

2. New annex type (obicompactvec, alongside bitmatrix.rs): a 4-bit- per-slot packed array (the presence mask above), one per partition — same on-disk shape family as PersistentBitMatrix's Packed variant, but simpler (no per-genome columns, a single derived read-only value per slot). As actually implemented (obicompactvec::siblingannex): not truly bit-packed — 1 byte/slot, 4 presence bits + the minorant bit (see the erratum above) in the low 5 bits, 3 unused. Deliberately simpler for a first implementation; packing to 5 bits/slot is a pure storage-density follow-up, not a behavioural change, still not done as of this note.

   Superseded 3-bit design (historical)

   3 bits, storing minorant status alongside sibling count directly, since it came for free from the same lookups (point 3 below) — 5 real states (not-minorant; minorant with 0/1/2/3 siblings) fit in 3 bits (8 states, 3 unused). This let the future SNP sweep discard a non-minorant slot instantly, with no lookup at all. The otherwise-unreachable combination "not-minorant + 0 siblings" (0 siblings always implies minorant) doubled as the "not yet computed" sentinel. Replaced by the 4-bit mask above, which subsumes this benefit (minorant still derivable, now for free at read time rather than stored) while also fixing the "which variant" blindness.

3. Computation pass (obikindex, new siblings.rs): one obipipeline run per layer, iterated sequentially over the index's layers — settled after two false starts, worth recording both.

   • False start 1: "fully parallel over every partition/slot at once, no ordering at all". Correctness is fine with this (sibling count and minorant are order-independent, unlike the old q >= p dedup they replace), but it reintroduces, at a larger scale, exactly the memory-blowup the original Step 2 sweep's sequential-outer-loop constraint existed to prevent: scattering every source partition at once multiplies the in-flight outgoing-query volume by the number of partitions.
   • False start 2: push the layer loop itself into the pipeline (source = the index's layers, a first Flat stage expands each layer into its k-mers). obipipeline's scheduler already bounds memory on its own — it dispatches every item through a shared worker pool at each stage boundary (scheduler.rs:217-372, dispatch() into a common worker_tx queue, any free worker picks up any pending item; not "one worker owns a chunk end to end"), with a biased Select that prioritises draining items already advanced in the chain over admitting new source items (scheduler.rs:271-282: stage results outrank the source, "vider le pipeline en priorité" / "dernier recours" for new data) — so bounded channel capacity plus this drain-first bias already caps in-flight work without any external sequential discipline. Correct, but it means k-mers from several layers can be completing concurrently, so the sink would need to track several open per-layer annex-file writers at once — real, avoidable complexity.
   • Settled design: keep the layer loop external and sequential — not for memory (the pipeline's own capacity/priority mechanism already provides that, for free, regardless), but so each pipeline run's sink targets exactly one layer's annex file, no concurrent multi-writer bookkeeping. Per layer: source = that layer's distinct k-mers; a Flat (1->N) stage generates the 3 central variants of a k-mer, each tagged with its origin (local slot); a transform stage routes each variant to its target partition (unchanged per-k-mer minimiser); a transform stage performs the lookup (existence-only — find_slot hit/miss, cheaper than the SNP sweep's full column fetch); a final stage/sink folds each answer into its origin's running state (below) and, once a layer's k-mers are all resolved, flushes the completed array to that layer's annex file. Many small, single-purpose stages on purpose, to let the scheduler interleave them finely across many in-flight items — this deliberately does not mirror how obipipeline is used elsewhere today: query.rs's process_chunk lumps parse+route+query+serialise into one closure (query.rs:325,743-758), and scatter.rs only pipelines file- reading/superkmer construction, routing partitions afterwards in a plain sequential loop (KmerPartition::write_batch, partition.rs:140) — both under-use the fine-grained scheduling the mechanism offers, so they are not precedents to copy, only existing (and arguably improvable, out of scope here) usages. Cross-partition lookups (querying another layer's MPHF for a variant) remain necessary as before — only the output side is kept single-layer.
   • Reconciliation: processed at the granularity of one answer batch per destination partition, not one source k-mer at a time — this is a proper shuffle, not a per-k-mer wait. Each source partition p holds a small array of running states (minorant = true, siblings = 0), one per local slot, initialised at scatter time and persisting across however many destination-partition batches answer it (up to 3, one per variant, not necessarily all from the same q). Every scattered query carries an origin tag (source partition + local slot) so its answer can be routed back. When target partition q returns its batch (all answers for every query that named q, regardless of which source k-mer or which source partition they came from), that batch is walked once, locally, and each answer updates — via its origin tag — the matching entry in its source partition's array: a miss changes nothing; a hit does siblings += 1, and if the found sibling's own encoding is smaller than the source's, minorant = false. A given source k-mer's state is final only once every destination batch concerning it has been folded in; its partition's array is flushed to the persistent annex once complete. Commutative per entry, so the order in which destination batches arrive and get folded in doesn't matter.

   Open optimisation, not adopted yet — real tradeoff, not a free win. Since looking up sibling y from x's visit already yields everything needed to fill y's own annex entry too, one visit per family could in principle replace one visit per observed family member — cutting this pass's cost roughly by the average family size instead of paying 3 * N_distinct regardless. But it means threads processing different source k-mers can end up writing the same sibling's slot concurrently — the fully independent, ownership-free parallelism of the plan above is deliberately traded away for this gain. It stays safe only because the computed value for a given slot is deterministic regardless of who computes it, so redundant concurrent writes converge to the same value — correct as long as each write is atomic, no locking needed — but it is a real design complexity increase over "every member redoes its own 3 lookups independently," not a strict improvement to adopt by default.

4. Trigger and caching (obikindex::KmerIndex/distance.rs): compute lazily on first phylo call for an SNP-family metric against a given index; check for an existing annex file first (mirrors PersistentBitMatrix::open()'s auto-detect-and-fall-back, bitmatrix.rs:264-287); if absent, run step 3 and persist; if present, mmap and reuse.

5. Invalidation. merge and filter always write to a fresh --output directory (obikindex/src/merge.rs, obikmer/src/cmd/filter.rs) so a re-merged/re-filtered index simply has no annex yet — nothing to invalidate. select --in-place (select_layer.rs:139-235) mutates columns of an existing index without changing its location, which can change sibling counts without removing rows — it must explicitly delete any cached annex for that index as part of its in-place rewrite.

6. Testing: hand-built tiny indexes with known sibling counts (0-3); order-independence (recompute twice on a static index, identical result, given the fully-parallel no-ownership design); invalidation (annex absent/correctly recomputed after select --in-place); once Step 2's sweep is fixed, a regression check that sibling_count == 0 slots are never looked up cross-partition during the sweep.

Cost: 3 * N_distinct existence-only lookups, computed once per index state and amortised over every subsequent phylo call that reuses the cached annex — cheaper per-lookup than the sweep itself (hit/miss only, no column fetch).

Step 3 — finalisation (obikindex)

From the global SnpTally alone (diagonal and off-diagonal both populated by the sweep, see Step 1/2 — no dependency on the external shared_kmers matrix), derive n x n distance matrices, each a pure function of the accumulated counts (same shape as jaccard_to_mash):

• p_hat[i,j] = SNP / (SNP + shared)
• Jukes-Cantor, Kimura-2P (from P, Q), optionally LogDet (needs the diagonal + base-composition margins).

Guard the singularities (p >= 3/4 for JC, 1-2P-Q <= 0 or 1-2Q <= 0 for K2P) by clamping to a max distance, as jaccard_to_mash clamps J <= 0.

Step 4 — surfacing (obikindex + obikmer CLI)

These metrics do not fit DistanceMetric's current LayeredStore-partial dispatch (they need the cross-partition sweep and produce a different intermediate). Two options, to decide:

• (a) New DistanceMetric variants (Pdistance, JukesCantor, Kimura2P, LogDet) whose KmerIndex::distance arm calls the sweep (snp.rs) instead of the partial path, still returning DistanceOutput. Keeps one CLI surface (--metric jukes-cantor), at the cost of a branch in distance() that ignores the LayeredStore it built.
• (b) A dedicated pathway (KmerIndex::snp_distance) and a distinct CLI entry, if mixing a cross-partition sweep into the partition-local phylo command is judged architecturally muddy.

Recommendation: (a) for user ergonomics (all pairwise distances under phylo, all feeding NJ/UPGMA/--shared-kmers unchanged), but compute the sweep lazily only when an SNP-family metric is requested, so the existing metrics keep their partition-local fast path untouched.

Step 5 — subsampling flag

Add --snp-sample <fraction> (or a bottom-s hash threshold): restrict the source-k-mer enumeration in Step 2 to seq_hash(kmer) < threshold. Divides lookups proportionally; p_hat is unbiased. Off by default (exact).

Testing

• Primitive unit tests: central_canonical_neighbors on hand-checked k-mers incl. palindrome-boundary cases; lone-k-mer minimizer against RollingStat's incremental result on the same k-mer.
• End-to-end tiny index: two 1-genome indexes differing by a handful of known isolated SNPs (transitions and transversions placed by hand), assert exact SNP, P, Q counts and the resulting JC/K2P values.
• Dedup invariant: assert the tally is identical regardless of genome/ partition order and that no pair is double-counted (compare against a brute-force all-pairs reference on a small index).
• Subsampling: p_hat within sampling error of the exact run.

Suggested phasing

1. Step 0 primitives + their unit tests (self-contained, no distance wiring). This also unblocks the long-declared-but-unimplemented query --mismatch (obikmer/src/cmd/query.rs:676, currently a warning), which needs the same neighbour + routing machinery.
2. SnpTally + finalisation math with a brute-force (non-swept) reference backend, validated on a tiny index.
3. The real per-partition sweep (Step 2) behind the same finalisation; assert it matches the brute-force backend.
4. CLI surfacing (Step 4a) and NJ/UPGMA integration (already generic over the matrix).
5. Subsampling (Step 5).

Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion)

Distinct from the SNP/Sankoff calibration above (which conditions on conserved-flank comparability), this addresses the coarser presence/absence signal of a family itself, for genome-skim/reduced-representation inputs where per-genome coverage is often <1x, sometimes <0.5x. At that depth, non-detection of a locus is dominated by Poisson sampling failure, not by true loss — even for a single-copy locus, P(non-detection | present) is non-negligible. Treating "absent" as a real 0-state therefore risks clustering genomes by shared undersampling rather than shared ancestry.

Parsimony. Loss and gain are the two directions of the same transition, so they cannot be freed asymmetrically without a biological justification for trusting one direction and not the other — none exists here. Freeing both directions to zero cost makes the character contribute nothing to any tree's score regardless of topology, which is operationally equivalent to just coding the state as ambiguous (?) rather than 0/1. ? is the more honest encoding when downstream ancestral-state reconstruction is wanted: cost-zero Fitch/Sankoff still assigns internal states arbitrarily among equal-cost options, whereas ? makes the "we don't know" explicit.

ML. The correct analogue is not a tuned asymmetric loss/gain rate but literal missing-data coding at the tip: an ambiguous state (all-states tip likelihood vector, e.g. (1,1) for a 2-state character) so Felsenstein pruning marginalises over it and the character contributes likelihood 1 (neutral) for that taxon — the same mechanism already used for the N IUPAC-ambiguous case in the 16-state Sankoff encoding (see "IQ-TREE custom-model format" above). No rate parameter to justify or estimate.

Open item, not yet decided. No per-locus/per-genome confidence distinction is implemented: the current plan is to code every non-detected locus as missing/ambiguous, uniformly, without a depth-based threshold for "confidently absent." Deferred because at typical skim coverage the confidently-absent case is rare enough that a uniform treatment is defensible as a first pass. Known consequence to watch for: a genome with very low overall coverage ends up with a mostly-? row and constrains the tree weakly — will likely need a minimum-detected-loci-per-genome filter upstream (exclude or flag under-covered genomes) before this becomes a correctness problem rather than a power problem. Not designed yet.

--free-loss implemented (2026-08-13)

obikmer phylo --sankoff --free-loss (with --tnt/--phyg/--iqtree) recodes ∅ (a family's non-detection in a genome) to ? instead of the default ordinary, costed 16th state — the resolution of the discussion above. ?, not -: - still carries gap/indel semantics in TNT/PhyG/ IQ-TREE (see "∅ stays an ordinary 16th state" above, and the reasoning that motivated it in the first place — the RAxML-era failure where - was silently swallowed as missing data was undesired back when ∅ was meant to be fully informative; --free-loss is the first case where that same tool behaviour is exactly what's wanted). ? is the standard missing/ambiguous symbol these tools already recognise without needing to be declared in the custom alphabet, so parsimony treats it as "compatible with any state" (Sankoff/Fitch's own missing-data extension — no per-transition cost edit needed) and Felsenstein pruning marginalises it for IQ-TREE's ML.

Implementation is a recoding at each export's existing --normalisation boundary, not a change to matrix calibration: write_sankoff_alignment_fasta (sankoff.rs, feeds --phyg too via its prefasta:), write_sankoff_tnt (tnt.rs, ? written directly, bypassing TNT_STATE_SYMBOL — there is no state index for it), and compact_alphabet/write_iqtree_alignment (iqtree.rs, ∅ excluded from the compact alphabet entirely rather than assigned a symbol). The calibrated cost matrix itself is unchanged — its ∅-row/column simply becomes unused once no sequence contains that state.

Unverified, flagged rather than assumed (this project's own standard for tool-specific claims — see the TNT/PhyG/IQ-TREE syntax facts above, all checked against the real binaries, not the manuals alone):

• PhyG auto-adds its own (n+1)-th gap/indel row to the tcm: matrix (phyg.rs); with --free-loss the reused _sankoff.fasta now contains real ? characters PhyG did not have to interpret before. Not confirmed against the local PhyG binary that ? is read as "any state" rather than routed through that gap row, or otherwise mishandled.
• IQ-TREE's exact ?/missing-data semantics for --seqtype MORPH with a custom-file model plus +ASC — assumed to marginalise via Felsenstein pruning by the general convention, not yet re-verified against the local iqtree3 build the way the rest of the custom-model mechanism was.

Both should be checked against a real run (small index) before --free-loss results are trusted for TNT/PhyG/IQ-TREE, the same way every other tool-specific behaviour in this file was established.

TNT confirmed (2026-08-13); IQ-TREE +ASC constant-column fix (2026-08-13)

TNT. Real run, 20-genome benchmark index, --free-loss --tnt: ? present for all 20 taxa in both sankoff.fasta and sankoff.tnt, loaded with no syntax error and no triangle-inequality warning, mult search completed, sankoff.tre holds 20 equally-parsimonious trees. Topology is exactly what the theory predicts: strain-level clades unaffected (the 4 E. coli, the Klebsiella trio, the Salmonella pairs — driven by real substitution signal, not gene-content) but resolution among the deeply divergent lineages (Candidozyma, Saccharolobus, Wolbachia, Acidobacterium, Opitutus, Bacillus/Shouchella) is weak and unstable across the 20 tied trees — the expected trade-off of discarding the gene-content/dropout signal that used to help place exactly those lineages. TNT's own handling of ? — read as "compatible with any state" — no longer a flagged assumption.

IQ-TREE: +ASC broke on real data — non-informative columns survive --free-loss, distinct bug from the --exclude-genome one above. snp_pseudo_alignment's "variable family" test (family_size() >= 2) is a whole-annex property, blind to any single column's actual calls. Once non-detection is hidden as ? (ignored for constancy) rather than counted as its own real state, a family that is genuinely variable somewhere in the dataset can still have only one distinct call among the genomes that actually carry it in one particular column, with the rest ? — constant once missing data is excluded, which is exactly what +ASC forbids. Same underlying failure mode as --exclude-genome's (documented above, "Two consistency bugs found and fixed post-implementation"), same remedy: rescan columns and drop the ones that are constant once - is ignored — but a different trigger (hiding cells vs. dropping whole rows), so a separate fix was needed rather than reusing drop_excluded directly. Implemented as iqtree::drop_ascertainment_noninformative, run only under --free-loss, only on IQ-TREE's own copy of the alignment (--tnt/--phyg have no no-invariant-site requirement and must not lose those columns). Not yet re-run against the local iqtree3 binary to confirm +ASC now accepts the filtered alignment — next verification step, alongside the still-open PhyG ?-vs-gap-row question above.

--family-overlap/--min-shared-family, and a second free_loss bug: cardinality-transition cost not zeroed (2026-08-13)

Diagnostic added: --family-overlap. Real per-genome-pair count of shared variable families (columns where neither genome is ∅), diagonal kept (degenerates to "genome's own total family count", not a separate computation — see obikmer/src/cmd/phylo/family_overlap.rs). On the 20-genome benchmark this predicted the --tnt/--iqtree rogue-placement pattern exactly: Candidozyma shares 0 families with all 19 others (hence its ~0-length branch, placed anywhere); the mean off-diagonal shared-family count per genome splits cleanly below/above ~1000 (Candidozyma 0, Saccharolobus 0.9, Opitutus 20.5, Wolbachia 37.1, Acidobacterium 41.2, Shouchella 49.4, Bacillus 54.5, Proteus 629, Yersinia 907.5 — all below 1000 — vs. Klebsiella ~5000-6000, Salmonella ~8000-17000, E. coli ~19000-24000, all comfortably above). --min-shared-family N automates exclusion below that threshold, reusing the --exclude-genome machinery. Restricting to the 3 genera above threshold (11 genomes: 4 E. coli, 3 Klebsiella, 4 Salmonella) gives 3 cleanly separated monophyletic clades in both TNT and IQ-TREE (without +ASC, see below) — confirms the family-overlap statistic, not just the tree topology, as the right lens for this failure mode. N=1000 is benchmark-specific (huge cross-domain divergence), not a universal constant — recompute --family-overlap and read its own gap on any other dataset.

Second free_loss bug, found by re-deriving the design rather than by a failing run: the cardinality-transition cost wasn't zeroed. --free-loss neutralises ∅ (whole-family non-detection) by recoding it to ?, so transitions to/from the empty state are never scored. But the calibrated 16-state cost matrix (obikindex::pairwise_cost_matrix) prices every cardinality change — not just X → ∅ — via a shared P_cardinality(|A|→|B|) factor (the same constant this project has called c_ctx throughout: "used everywhere a member is gained or lost", see "A concrete Sankoff cost matrix" above). So a transition between two non-empty states that differ in how many family members are detected (e.g. {A} -> {A,C}, gaining a second sibling) was still carrying the full calibrated gain/loss penalty under --free-loss, even though the same low/incomplete-coverage argument for not trusting whole-family absence applies identically one level down: whether a genome shows 1 vs 2 detected members of a family it does carry is exactly as vulnerable to sampling failure as whether the family was detected at all. Fixed by adding a free_loss: bool parameter to pairwise_cost_matrix that drops the P_cardinality factor from the log-probability sum entirely when set, leaving cost driven only by composition matching (shared-base retention, paired substitutions) — never by cardinality difference alone. Unit-tested (free_loss_ignores_cardinality_transition_probs): two wildly different cardinality models produce an identical cost matrix once free_loss is true. Not yet re-run against a real index to see how much this changes the resulting tree beyond the unit-level guarantee.

+ASC misbehaves on closely-related taxa, separately from either bug above. Restricting to the 3-genus, 11-genome subset and running --iqtree with +ASC converges suspiciously fast compared to the same run without +ASC (which takes the expected ~100 NNI iterations, ~1m38s, converges to a stable score, total tree length 7.147, and gives the clean 3-clade topology). Not yet root-caused with a real +ASC log in hand (the one run overwrote its own log under the shared --prefix, see the operational note below) — working hypothesis, not confirmed: +ASC's correction divides by P(site variable | current model, current branch lengths), recomputed every iteration; on closely-related strains (short branches), a plausible model predicts most sites should be constant, so P(variable) can be small enough that the correction becomes numerically unstable. Needs a dedicated run (own --prefix, log kept) before this is more than a hypothesis.

Operational note, not a code bug: IQ-TREE's --undo/checkpoint reuse across different alignments crashes. Reusing iqtree.ckp.gz (and the tree it stores) from a run on one taxon set against a regenerated iqtree.fasta with a different taxon count crashes with Assertion 'i >= 0 && i < (int)seq_names.size()' failed in Alignment::getSeqName / PhyloTree::assignLeafNames — --undo tries to graft a stored tree with the old leaf count onto the new, smaller alignment. Not an obikmer bug; avoid by using a fresh --prefix (or -redo) whenever the underlying alignment's taxon set changes, never --undo across them.

_iqtree_states.csv: compact-symbol traceability (2026-08-15)

User-reported suspicion: --iqtree --free-loss state frequencies "don't sum to 1". Verified both by code trace and unit test (obikmer/src/cmd/phylo/iqtree.rs, free_loss_excludes_absent_state_and_freq_sums_to_one) — not a bug: compact_alphabet's counting loop continues on free_loss && b == b'-' before recoding/counting, so the absent state never enters occurs/counts, and freq[i] = counts[old]/total sums to 1.0 by construction (total is defined as the sum over exactly the same states). Confirmed against real data the user provided: a pasted _iqtree.model frequency line summed to 1.000001 — a 6-decimal print-rounding artifact (format!("{p:.6}") × 15 values), not a computation error; IQ-TREE's own 3-decimal-rounded report of the same values summed to exactly 1.000.

Investigating that report surfaced a real, separate gap: nothing mapped _iqtree.model/_iqtree.fasta's compact state symbols (0-9A-F, IQ-TREE renumbers away unused states from the full 16, see this module's own "--iqtree... state count" discussion) back to the canonical 16-state alphabet (STATE_SYMBOL, sankoff.rs) — so a pattern like "compact state 0 has zero exchangeability with every other state" (R(a,b) = exp(-cost(a,b)) = 0 for an entire row/column) couldn't be traced back to which real state that was, or whether it was expected (a state combination that was simply never observed alongside anything else in the calibration data, giving cost = -ln(0) = +∞) or a genuine problem.

Fix: write_iqtree_states_csv (iqtree.rs) writes <prefix>_iqtree_states.csv — one row per surviving state, iqtree_symbol,canonical_symbol,frequency, frequency at full f64 precision (not the model file's truncated 6 decimals). Written alongside _iqtree.model/_iqtree.fasta from the same CompactAlphabet both already use, so there is no risk of the three files disagreeing. Covered by states_csv_maps_compact_symbols_back_to_canonical_ones.

Root cause of the zero-exchangeability pattern found and fixed (2026-08-15): premature exp() in pairwise_cost_matrix underflowed merely-tiny probabilities to exactly 0.0. The user also reported iqtree3 emitting "Numerical underflow for lh-derivative" warnings on the same run — a real signal, traced to obikphylo/src/cardcomp.rs's pairwise_cost_matrix, not to the frequency computation (which is a plain, safe f64 division, never close to underflow at any realistic scale). The function already accumulated log_p in log-space (correct), but then row-normalised by exponentiating each cell first (raw[a][b] = log_p.exp()) and summing the results — f64::exp hard underflows to exactly 0.0 for any input below roughly -709, which a sum of several individually-small-but-nonzero probability factors (composition/cardinality terms, best_pairing_cost's pairing terms) can reach easily on real, skewed calibration data. Once raw[a][b] was exactly 0.0, normalisation and -ln turned a merely tiny probability into a +∞ cost indistinguishable from a literally unobserved transition (p == 0.0 exactly, e.g. p_comp[i][j] never once tallied) — conflating two different things: "never observed" (should be +∞, a correct MLE result) and "observed, but the joint probability of this multi-step transition is extremely small" (should be a large finite cost).

Fix: row-normalise via log-sum-exp instead of exponentiating first — row_max = max_b(log_p[a][b]), log_sum = row_max + ln(Σ_b exp(log_p[a][b] - row_max)) (every shifted term is in (0,1], so this never underflows for a finite log_p[a][b]), then cost[a][b] = log_sum - log_p[a][b] directly — no intermediate probability is ever materialised. This falls out of IEEE 754 arithmetic without a special case: a genuinely-unobserved factor (log_p[a][b] == -∞, from the existing if p > 0.0 {...} else { NEG_INFINITY } guards already in the log-accumulation loop) still yields cost = +∞ exactly (finite − (−∞) = +∞), preserving the correct semantics for that case, while every merely-tiny-but-nonzero transition now gets a large but finite cost. A degenerate all--∞ row (a state with literally zero probability of transitioning to anything, row_max == -∞) is guarded explicitly to avoid a -∞ − (-∞) = NaN in the log-sum-exp itself. Covered by cardcomp::tests::underflow_prone_transition_gets_finite_cost_not_infinite (all off-diagonal composition probabilities set to 1e-200, well past where the old exp()-first code would have underflowed to 0.0, cost asserted finite). Every pre-existing cardcomp test still passes unchanged (numerically identical results when no underflow occurs).

References

The Mash mutation-rate model this discussion contrasts with: [@Mash-distances-doc; @Fan2015-mash-formula].