From 96542018851aeff4ef7981035c825e5ec91a6dd3 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 15 Aug 2026 21:30:54 +0200 Subject: [PATCH] feat: introduce _iqtree_states.csv for compact symbol mapping Generates a new CSV output that maps IQ-TREE's compact state symbols to canonical states alongside full-precision empirical frequencies. Updates documentation to clarify that state frequencies sum to 1.0 by design and documents conditional behavior under `--free-loss`. Includes unit tests verifying absent state exclusion, frequency summation, and CSV structure. Also restricts entropy annex resolution to non-monomorphic minorants to eliminate redundant per-genome checks. --- DevDoc/architecture/siblings/index.html | 21 +++- .../theory/evolutionary_distances/index.html | 58 ++++++++++ DevDocMD/theory/evolutionary_distances.md | 41 +++++++ UserDocMD/usage/phylo.md | 3 + doc/usage/phylo/index.html | 7 ++ src/obikmer/src/cmd/phylo/iqtree.rs | 104 +++++++++++++++++- 6 files changed, 230 insertions(+), 4 deletions(-) diff --git a/DevDoc/architecture/siblings/index.html b/DevDoc/architecture/siblings/index.html index 816ff532..8062296e 100644 --- a/DevDoc/architecture/siblings/index.html +++ b/DevDoc/architecture/siblings/index.html @@ -2051,13 +2051,28 @@ but indexed by family_idx (every minorant of the layer, monomorphic included — the same numbering Selection/scan_layer_families already use), one f32 entropy15 value per entry, -1.0 sentinel for monomorphic/ not-yet-computed. First use of --entropy/--entropy-sd on an index -pays a one-time cost (ensure_entropy_annexes in entropy.rs: a full, -unsampled Selection::All scan, resolving every non-monomorphic -minorant's genome_mask once to compute and persist its entropy) — every +pays a one-time cost (ensure_entropy_annexes in entropy.rs) — every later run (any μ/σ, any command) reads the file positionally, no re-scan, restoring the usual Selection::Some "skip resolving excluded families" speedup that a naive "weigh during the resolving scan" design would have permanently forfeited.

+

Bug found and fixed (2026-08-15): ensure_entropy_annexes scanned with +Selection::All instead of bounding to non-monomorphic minorants. +Monomorphism (family_size() < 2) is knowable directly from the annex +bits alone, no per-genome resolution needed — but the original +implementation called scan_layer_families with Selection::All +anyway, so fill_sub_matrix_carries (the expensive per-genome +resolution) ran for every minorant, ~98% of which are monomorphic +(measured elsewhere in this doc) and had their genome_mask immediately +discarded once the callback checked family_size() < 2. Fixed by adding +[subsample::non_monomorphic_selection_layer] — a cheap, annex-only, +non-sampling pass (same shape as reservoir_sample_layer, but keeping +every non-monomorphic minorant's family_idx instead of a bounded +reservoir) — and passing Selection::Some(&eligible) instead of +Selection::All, so the expensive resolution now runs only for the ~2% +of minorants that can actually produce a real entropy value. A +debug_assert!(mask.family_size() >= 2, ...) inside the +scan_layer_families callback guards the invariant.

Resolved: the existing hard "non-monomorphic minorant" eligibility filter stays a hard gate upstream of the Gaussian weighting — only qualifying families ever get a stored entropy value or a weighted draw.

diff --git a/DevDoc/theory/evolutionary_distances/index.html b/DevDoc/theory/evolutionary_distances/index.html index c3b758ce..c56c4b9e 100644 --- a/DevDoc/theory/evolutionary_distances/index.html +++ b/DevDoc/theory/evolutionary_distances/index.html @@ -887,6 +887,17 @@ + + +
  • + + + + _iqtree_states.csv: compact-symbol traceability (2026-08-15) + + + +
  • @@ -2095,6 +2106,17 @@ +
    • + +
  • + + + + _iqtree_states.csv: compact-symbol traceability (2026-08-15) + + + +
  • @@ -4048,6 +4070,42 @@ with a different taxon count crashes with Assertion 'i >= 0 && 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.

    +

    The zero-exchangeability pattern itself (state 0 in the user's report, +frequency 3.26%, R=0 with every other state) is not yet explained — +plausibly a genuinely unobserved transition in the calibration +(cardinality_transitions/composition_transitions count 0 for every +pair involving it), which is a legitimate, if numerically extreme, result +of -ln(0), not necessarily a bug — not investigated further.

    References

    The Mash mutation-rate model this discussion contrasts with: (Fan et al. 2015; Marbl Lab 2026)1 2.

    diff --git a/DevDocMD/theory/evolutionary_distances.md b/DevDocMD/theory/evolutionary_distances.md index 9405adfc..23d1805d 100644 --- a/DevDocMD/theory/evolutionary_distances.md +++ b/DevDocMD/theory/evolutionary_distances.md @@ -2056,6 +2056,47 @@ 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 `continue`s 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 +`_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`. + +The zero-exchangeability pattern itself (state 0 in the user's report, +frequency 3.26%, `R=0` with every other state) is not yet explained — +plausibly a genuinely unobserved transition in the calibration +(`cardinality_transitions`/`composition_transitions` count `0` for every +pair involving it), which is a legitimate, if numerically extreme, result +of `-ln(0)`, not necessarily a bug — not investigated further. + ## References The Mash mutation-rate model this discussion contrasts with: diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md index 2d068200..8d0fe5a5 100644 --- a/UserDocMD/usage/phylo.md +++ b/UserDocMD/usage/phylo.md @@ -210,6 +210,7 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below): | `_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search | | `_iqtree.model` | `--iqtree` | IQ-TREE model file | custom ML substitution model | | `_iqtree.fasta` | `--iqtree` | FASTA | alignment recoded for that model | +| `_iqtree_states.csv` | `--iqtree` | CSV table | maps `_iqtree.model`/`_iqtree.fasta`'s compact state symbols back to `_sankoff_matrix.csv`'s alphabet | **`_sankoff_matrix.csv`** — header `state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N`: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), `0` standing for the empty/absent state (not `-`, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost $-\ln P(a,b)$, 4 decimals. @@ -235,3 +236,5 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below): ``` iqtree3 -s _iqtree.fasta --seqtype MORPH -m _iqtree.model+ASC --prefix _iqtree -T AUTO ``` + +**`_iqtree_states.csv`** (`--iqtree`) — one row per state actually kept in `_iqtree.model`/`_iqtree.fasta` (header `iqtree_symbol,canonical_symbol,frequency`): `iqtree_symbol` is the compact `0-9A-F` symbol as written in those two files, `canonical_symbol` is the matching `_sankoff_matrix.csv` state, `frequency` is that state's empirical frequency at full precision (`_iqtree.model`'s own frequency line is rounded to 6 decimals). Under `--free-loss`, absent (`0`/`?`) is never a kept state, so it never appears here. Use this file to identify which real state a given row/column of `_iqtree.model`'s matrix corresponds to — e.g. to check whether a state showing zero exchangeability with everything else is expected (a state combination that never co-occurs with anything else in this data) or worth investigating further. diff --git a/doc/usage/phylo/index.html b/doc/usage/phylo/index.html index 3a47990f..e19e2aca 100644 --- a/doc/usage/phylo/index.html +++ b/doc/usage/phylo/index.html @@ -2014,6 +2014,12 @@ FASTA alignment recoded for that model + +<prefix>_iqtree_states.csv +--iqtree +CSV table +maps _iqtree.model/_iqtree.fasta's compact state symbols back to _sankoff_matrix.csv's alphabet +

    _sankoff_matrix.csv — header state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), 0 standing for the empty/absent state (not -, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost \(-\ln P(a,b)\), 4 decimals.

    @@ -2048,6 +2054,7 @@

    _iqtree.fasta (--iqtree) — alignment recoded to that same compact 0..k-1 alphabet (symbols 0-9A-F). Under --free-loss, non-detection becomes ? and columns left non-informative once missing calls are ignored are dropped first (required for +ASC). Run with: 

    iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO

    +

    _iqtree_states.csv (--iqtree) — one row per state actually kept in _iqtree.model/_iqtree.fasta (header iqtree_symbol,canonical_symbol,frequency): iqtree_symbol is the compact 0-9A-F symbol as written in those two files, canonical_symbol is the matching _sankoff_matrix.csv state, frequency is that state's empirical frequency at full precision (_iqtree.model's own frequency line is rounded to 6 decimals). Under --free-loss, absent (0/?) is never a kept state, so it never appears here. Use this file to identify which real state a given row/column of _iqtree.model's matrix corresponds to — e.g. to check whether a state showing zero exchangeability with everything else is expected (a state combination that never co-occurs with anything else in this data) or worth investigating further.

    diff --git a/src/obikmer/src/cmd/phylo/iqtree.rs b/src/obikmer/src/cmd/phylo/iqtree.rs index 6d659630..b774ce59 100644 --- a/src/obikmer/src/cmd/phylo/iqtree.rs +++ b/src/obikmer/src/cmd/phylo/iqtree.rs @@ -5,7 +5,7 @@ use obifastwrite::{JsonVal, write_record}; use obikphylo::siblings::SnpAlignment; use tracing::info; -use super::sankoff::state_index_table; +use super::sankoff::{state_index_table, STATE_SYMBOL}; // ── Sankoff-calibrated data → IQ-TREE custom ML model + recoded alignment ── // @@ -145,6 +145,30 @@ fn compact_alphabet(alignment: &SnpAlignment, free_loss: bool) -> CompactAlphabe CompactAlphabet { old_to_compact, compact_to_old, freq } } +/// Write `_iqtree_states.csv`: the mapping from IQ-TREE's own +/// compact state symbols (`0..9A-F`, what actually appears in +/// `_iqtree.fasta`/`_iqtree.model`) back to the canonical 16-state +/// alphabet (`STATE_SYMBOL` — the same one `_sankoff_matrix.csv` is +/// indexed by), plus each state's empirical frequency at full precision +/// (`_iqtree.model`'s own frequency line is truncated to 6 decimals). +/// Without this file, a compact index in `_iqtree.model`'s `R`/`π` output +/// (e.g. "state 0 has zero exchangeability with everything else") can't be +/// traced back to which real state that is. +fn write_iqtree_states_csv(alphabet: &CompactAlphabet, output: &Option) -> String { + let path = output.as_ref() + .map(|p| format!("{}_iqtree_states.csv", p.display())) + .unwrap_or_else(|| "iqtree_states.csv".into()); + let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| { + eprintln!("error creating {path}: {e}"); + std::process::exit(1); + })); + writeln!(f, "iqtree_symbol,canonical_symbol,frequency").unwrap(); + for (compact, &old) in alphabet.compact_to_old.iter().enumerate() { + writeln!(f, "{},{},{}", IQTREE_STATE_SYMBOL[compact], STATE_SYMBOL[old as usize], alphabet.freq[compact]).unwrap(); + } + path +} + /// Write the `R` (exchangeability) + `π` (frequencies) model file IQ-TREE's /// `-m +ASC` reads. Returns the path, so the caller can print a /// single combined "how to run this" message once the alignment is also @@ -234,6 +258,7 @@ pub(super) fn write_iqtree( }; let alphabet = compact_alphabet(alignment, free_loss); + let states_path = write_iqtree_states_csv(&alphabet, output); let model_path = write_iqtree_model(matrix, &alphabet, output); let (fasta_path, n_sites) = write_iqtree_alignment(alignment, labels, &alphabet, output, free_loss); @@ -243,8 +268,85 @@ pub(super) fn write_iqtree( .unwrap_or_else(|| "iqtree".into()); info!( "IQ-TREE alignment → {fasta_path} ({n_sites} sites, {} states)\n\ + IQ-TREE state mapping → {states_path}\n\ Run with:\n \ iqtree3 -s {fasta_path} --seqtype MORPH -m {model_path}+ASC --prefix {prefix_name} -T AUTO", alphabet.k() ); } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn free_loss_excludes_absent_state_and_freq_sums_to_one() { + // 3 genomes, 2 sites. Site 0: g1='A', g2='C', g3='-' (absent). + // Site 1: g1='-', g2='-', g3='G'. Under free_loss, every '-' must + // be excluded from the frequency count entirely (not folded into + // state 0) — reproduces a user-reported suspicion that state 0 + // ("absent") was still being counted despite being recoded to `?` + // (IQ-TREE's own missing symbol) in the alignment actually written. + let alignment = SnpAlignment { + sequences: vec![ + vec![b'A', b'-'], + vec![b'C', b'-'], + vec![b'-', b'G'], + ], + }; + + let alphabet = compact_alphabet(&alignment, true); + + assert!( + !alphabet.compact_to_old.contains(&0), + "state 0 (absent) must not appear in the compact alphabet under --free-loss, got {:?}", + alphabet.compact_to_old + ); + let sum: f64 = alphabet.freq.iter().sum(); + assert!((sum - 1.0).abs() < 1e-9, "frequencies must sum to 1, got {sum} ({:?})", alphabet.freq); + assert_eq!(alphabet.k(), 3, "A, C, G — 3 real states, `-` excluded"); + } + + #[test] + fn without_free_loss_absent_state_is_counted_normally() { + let alignment = SnpAlignment { + sequences: vec![ + vec![b'A', b'-'], + vec![b'C', b'-'], + vec![b'-', b'G'], + ], + }; + + let alphabet = compact_alphabet(&alignment, false); + + assert!( + alphabet.compact_to_old.contains(&0), + "state 0 (absent, recoded from '-') must be counted when --free-loss is off" + ); + let sum: f64 = alphabet.freq.iter().sum(); + assert!((sum - 1.0).abs() < 1e-9, "frequencies must sum to 1, got {sum} ({:?})", alphabet.freq); + } + + #[test] + fn states_csv_maps_compact_symbols_back_to_canonical_ones() { + // 'A' (state 1) and 'G' (state 4) occur, '-' (state 0) excluded by + // --free-loss — compact index 0 -> 'A', compact index 1 -> 'G'. + let alignment = SnpAlignment { + sequences: vec![ + vec![b'A', b'-'], + vec![b'-', b'G'], + ], + }; + let alphabet = compact_alphabet(&alignment, true); + let output = Some(std::env::temp_dir().join(format!("obikmer_test_iqtree_states_{}", std::process::id()))); + + let path = write_iqtree_states_csv(&alphabet, &output); + let csv = std::fs::read_to_string(&path).unwrap(); + std::fs::remove_file(&path).ok(); + let mut lines = csv.lines(); + assert_eq!(lines.next(), Some("iqtree_symbol,canonical_symbol,frequency")); + assert_eq!(lines.next(), Some("0,A,0.5")); + assert_eq!(lines.next(), Some("1,G,0.5")); + assert!(lines.next().is_none()); + } +}