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