From 128db645644266f6e2acd42014a1c3e5b5ea6bf6 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 15 Aug 2026 22:38:40 +0200 Subject: [PATCH] feat(phylo): add --iqtree-min-freq to filter rare nucleotide states Introduces --iqtree-min-freq (default 0.001) to treat low-frequency nucleotide states as missing data during IQ-TREE alignment generation when --free-loss is active. This triggers a recoding pass that folds rare states into the missing symbol, followed by non-informative site removal and alphabet recomputation to maintain output consistency. The change also adds Sankoff model configuration files and updates related tests and documentation. --- .../theory/evolutionary_distances/index.html | 65 +++++++ DevDocMD/theory/evolutionary_distances.md | 50 ++++++ UserDocMD/usage/phylo.md | 9 +- aloss_nofl_sankoff_params.yaml | 167 ++++++++++++++++++ aloss_sankoff_params.yaml | 167 ++++++++++++++++++ doc/usage/phylo/index.html | 33 +++- src/obikmer/src/cmd/phylo/args.rs | 19 ++ src/obikmer/src/cmd/phylo/iqtree.rs | 85 ++++++++- src/obikmer/src/cmd/phylo/mod.rs | 2 +- src/obikphylo/src/siblings/tests.rs | 166 ++++++++++++++++- 10 files changed, 756 insertions(+), 7 deletions(-) create mode 100644 aloss_nofl_sankoff_params.yaml create mode 100644 aloss_sankoff_params.yaml diff --git a/DevDoc/theory/evolutionary_distances/index.html b/DevDoc/theory/evolutionary_distances/index.html index 71d10aa8..e7627d12 100644 --- a/DevDoc/theory/evolutionary_distances/index.html +++ b/DevDoc/theory/evolutionary_distances/index.html @@ -898,6 +898,17 @@ + + +
  • + + + + --iqtree-min-freq: rare states still perturb IQ-TREE after the underflow fix (implemented 2026-08-15) + + + +
  • @@ -2117,6 +2128,17 @@ +
    • + +
  • + + + + --iqtree-min-freq: rare states still perturb IQ-TREE after the underflow fix (implemented 2026-08-15) + + + +
  • @@ -4140,6 +4162,49 @@ Covered by cardcomp::tests::underflow_prone_transition_gets_finite_cost_no 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).

    +

    --iqtree-min-freq: rare states still perturb IQ-TREE after the underflow fix (implemented 2026-08-15)

    +

    After the pairwise_cost_matrix fix above, the user re-checked +_iqtree_states.csv on the real phyloskims_sal_vac run: no more exact +zeros, but N (all 4 bases), V/H/D (3 bases each) sit at +8×10⁻⁵–4×10⁻⁴, and B (3 bases) just above 10⁻³ — all four +3-simultaneous-base states plus N, exactly the states the underflow fix +was papering over the symptom of. iqtree3 kept emitting "Numerical +underflow for lh-derivative" on this same run — a model with several +near-zero-frequency states is still numerically marginal for ML +optimisation even once no cost/exchangeability entry is exactly +infinite/zero.

    +

    Decided and implemented: extend --free-loss's existing "sampling +failure, not true signal" treatment of non-detection to rare +multi-base states too — biologically, a family showing 3–4 simultaneously +observed central bases across genomes is far more likely to be +assembly/detection noise than a genuinely ancestral, widely-preserved +polymorphism, the same reasoning already used for -.

    +

    Mechanism (obikmer/src/cmd/phylo/iqtree.rs, write_iqtree): a +second filtering pass, active only under --free-loss. After the first +compact_alphabet call (post ascertainment-filtering), any state with +freq < --iqtree-min-freq (default 0.001) has every occurrence recoded +to - (recode_symbols_as_absent, a new small helper) — reusing the +existing -→? machinery in compact_alphabet/write_iqtree_alignment +rather than duplicating it, since a state recoded to - flows through +exactly the same downstream path a genuinely-absent call already does. +drop_ascertainment_noninformative then re-runs on the recoded alignment +(a site can become non-informative once its only distinguishing call was +the now-recoded rare state), and compact_alphabet is recomputed a third +time on the result — this final alphabet is what _iqtree_states.csv/ +_iqtree.model/_iqtree.fasta all get built from, so the three files +stay consistent with each other by construction, same as before this +change.

    +

    Not a fixed-point loop (recoding one rare state can't, by construction, +lower another already-above-threshold state's frequency enough to also +cross the threshold, since removing calls only redistributes mass upward +among the survivors) — a single extra pass is sufficient.

    +

    --iqtree-min-freq 0 disables the filter (no state can have frequency +< 0), recovering the pre-this-change behavior exactly.

    +

    Covered by iqtree::tests::iqtree_min_freq_folds_rare_states_into_missing +(20 common A/C sites + one site introducing a single, otherwise-informative +M call at ~1/62 frequency, --iqtree-min-freq 0.05; asserts M absent +from the written _iqtree_states.csv and A/C still present). Full +workspace cargo test green.

    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 d2d227bf..d1cf4927 100644 --- a/DevDocMD/theory/evolutionary_distances.md +++ b/DevDocMD/theory/evolutionary_distances.md @@ -2132,6 +2132,56 @@ 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). +## `--iqtree-min-freq`: rare states still perturb IQ-TREE after the underflow fix (implemented 2026-08-15) + +After the `pairwise_cost_matrix` fix above, the user re-checked +`_iqtree_states.csv` on the real `phyloskims_sal_vac` run: no more exact +zeros, but `N` (all 4 bases), `V`/`H`/`D` (3 bases each) sit at +`8×10⁻⁵`–`4×10⁻⁴`, and `B` (3 bases) just above `10⁻³` — all four +3-simultaneous-base states plus `N`, exactly the states the underflow fix +was papering over the symptom of. `iqtree3` kept emitting "Numerical +underflow for lh-derivative" on this same run — a model with several +near-zero-frequency states is still numerically marginal for ML +optimisation even once no cost/exchangeability entry is *exactly* +infinite/zero. + +**Decided and implemented**: extend `--free-loss`'s existing "sampling +failure, not true signal" treatment of non-detection to rare +multi-base states too — biologically, a family showing 3–4 simultaneously +observed central bases across genomes is far more likely to be +assembly/detection noise than a genuinely ancestral, widely-preserved +polymorphism, the same reasoning already used for `-`. + +**Mechanism** (`obikmer/src/cmd/phylo/iqtree.rs`, `write_iqtree`): a +*second* filtering pass, active only under `--free-loss`. After the first +`compact_alphabet` call (post ascertainment-filtering), any state with +`freq < --iqtree-min-freq` (default `0.001`) has every occurrence recoded +to `-` (`recode_symbols_as_absent`, a new small helper) — reusing the +*existing* `-`→`?` machinery in `compact_alphabet`/`write_iqtree_alignment` +rather than duplicating it, since a state recoded to `-` flows through +exactly the same downstream path a genuinely-absent call already does. +`drop_ascertainment_noninformative` then re-runs on the recoded alignment +(a site can become non-informative once its only distinguishing call was +the now-recoded rare state), and `compact_alphabet` is recomputed a third +time on the result — this final alphabet is what `_iqtree_states.csv`/ +`_iqtree.model`/`_iqtree.fasta` all get built from, so the three files +stay consistent with each other by construction, same as before this +change. + +Not a fixed-point loop (recoding one rare state can't, by construction, +lower another *already-above-threshold* state's frequency enough to also +cross the threshold, since removing calls only redistributes mass upward +among the survivors) — a single extra pass is sufficient. + +`--iqtree-min-freq 0` disables the filter (no state can have frequency +`< 0`), recovering the pre-this-change behavior exactly. + +Covered by `iqtree::tests::iqtree_min_freq_folds_rare_states_into_missing` +(20 common A/C sites + one site introducing a single, otherwise-informative +`M` call at ~1/62 frequency, `--iqtree-min-freq 0.05`; asserts `M` absent +from the written `_iqtree_states.csv` and `A`/`C` still present). Full +workspace `cargo test` green. + ## References The Mash mutation-rate model this discussion contrasts with: diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md index 8d0fe5a5..7da90cc9 100644 --- a/UserDocMD/usage/phylo.md +++ b/UserDocMD/usage/phylo.md @@ -124,6 +124,7 @@ The first `phylo` run on a given index that uses `--entropy`/`--entropy-sd` pays | `--tnt` | off | Also write a TNT script (implies `--sankoff`) | | `--phyg` | off | Also write PhyG input files (implies `--sankoff`) | | `--iqtree` | off | Also write an IQ-TREE custom model and alignment (implies `--sankoff`) | +| `--iqtree-min-freq` | `0.001` | With `--iqtree --free-loss`: also treat as missing any state rarer than this in the alignment | | `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands) | ### The 16-state model @@ -232,9 +233,13 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below): **`_iqtree.model`** (`--iqtree`) — lower-triangular exchangeability matrix $R(a,b) = e^{-\text{cost}(a,b)}$ (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered `0..k-1`. -**`_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: +**`_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`); with `--iqtree-min-freq` also set (the default), any state rarer than that threshold is folded into the same `?` treatment, and non-informative columns are re-checked and dropped again after that. Run with: ``` 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. +**`_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 — nor does any state `--iqtree-min-freq` folded away for being too rare. 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. + +### Rare states and `--iqtree-min-freq` + +States that combine 3 or 4 central bases at once (IUPAC `V`/`H`/`K`.../`N`) are inherently rare — and, on real data, rare enough that they can make `iqtree3` itself numerically unstable ("Numerical underflow for lh-derivative" warnings, near-degenerate likelihood optimization). They're also more likely to be assembly/detection noise than genuine, widely-shared multi-way polymorphism, the same "sampling failure, not true signal" reasoning `--free-loss` already applies to non-detection. With `--free-loss` set, `--iqtree-min-freq` (default `0.001`, i.e. one in a thousand) extends that same missing-data treatment to any state below this frequency, not just absence. Check `_iqtree_states.csv` to see exactly which states survived and at what frequency; set `--iqtree-min-freq 0` to disable this and keep every state that occurs at all (the old behavior). Has no effect without `--free-loss` — there is no missing-data symbol to fold rare states into otherwise. diff --git a/aloss_nofl_sankoff_params.yaml b/aloss_nofl_sankoff_params.yaml new file mode 100644 index 00000000..1e488fab --- /dev/null +++ b/aloss_nofl_sankoff_params.yaml @@ -0,0 +1,167 @@ +ratio_ceiling: 0.5 +cardinality_transitions: +- from: 0 + to: 0 + count: 122014779 + probability: 0.870049812090934 +- from: 0 + to: 1 + count: 17966272 + probability: 0.12811195254940888 +- from: 0 + to: 2 + count: 233683 + probability: 0.001666321505518981 +- from: 0 + to: 3 + count: 24109 + probability: 0.00017191385413811494 +- from: 0 + to: 4 + count: 0 + probability: 0.0 +- from: 1 + to: 0 + count: 17966272 + probability: 0.967023782579572 +- from: 1 + to: 1 + count: 602798 + probability: 0.03244523972983381 +- from: 1 + to: 2 + count: 9562 + probability: 0.0005146688978673966 +- from: 1 + to: 3 + count: 303 + probability: 0.00001630879272681669 +- from: 1 + to: 4 + count: 0 + probability: 0.0 +- from: 2 + to: 0 + count: 233683 + probability: 0.9585500516842502 +- from: 2 + to: 1 + count: 9562 + probability: 0.039222603245442765 +- from: 2 + to: 2 + count: 438 + probability: 0.0017966429848885097 +- from: 2 + to: 3 + count: 105 + probability: 0.00043070208541847837 +- from: 2 + to: 4 + count: 0 + probability: 0.0 +- from: 3 + to: 0 + count: 24109 + probability: 0.9827171564831044 +- from: 3 + to: 1 + count: 303 + probability: 0.012350711286838137 +- from: 3 + to: 2 + count: 105 + probability: 0.004279949455834998 +- from: 3 + to: 3 + count: 16 + probability: 0.0006521827742224758 +- from: 3 + to: 4 + count: 0 + probability: 0.0 +- from: 4 + to: 0 + count: 0 + probability: 0.0 +- from: 4 + to: 1 + count: 0 + probability: 0.0 +- from: 4 + to: 2 + count: 0 + probability: 0.0 +- from: 4 + to: 3 + count: 0 + probability: 0.0 +- from: 4 + to: 4 + count: 0 + probability: 0.0 +composition_transitions: +- from: 'A' + to: 'A' + count: 0 + probability: 0.0 +- from: 'A' + to: 'C' + count: 0 + probability: 0.0 +- from: 'A' + to: 'G' + count: 0 + probability: 0.0 +- from: 'A' + to: 'T' + count: 0 + probability: 0.0 +- from: 'C' + to: 'A' + count: 0 + probability: 0.0 +- from: 'C' + to: 'C' + count: 96389 + probability: 0.9663832688335907 +- from: 'C' + to: 'G' + count: 1290 + probability: 0.012933368089671353 +- from: 'C' + to: 'T' + count: 2063 + probability: 0.020683363076737984 +- from: 'G' + to: 'A' + count: 0 + probability: 0.0 +- from: 'G' + to: 'C' + count: 1290 + probability: 0.005696948820201646 +- from: 'G' + to: 'G' + count: 222720 + probability: 0.9835848381669073 +- from: 'G' + to: 'T' + count: 2427 + probability: 0.010718213012891003 +- from: 'T' + to: 'A' + count: 0 + probability: 0.0 +- from: 'T' + to: 'C' + count: 2063 + probability: 0.007305266661709142 +- from: 'T' + to: 'G' + count: 2427 + probability: 0.008594223067362136 +- from: 'T' + to: 'T' + count: 277909 + probability: 0.9841005102709287 diff --git a/aloss_sankoff_params.yaml b/aloss_sankoff_params.yaml new file mode 100644 index 00000000..e30daac6 --- /dev/null +++ b/aloss_sankoff_params.yaml @@ -0,0 +1,167 @@ +ratio_ceiling: 0.5 +cardinality_transitions: +- from: 0 + to: 0 + count: 13019012 + probability: 0.8682962128540935 +- from: 0 + to: 1 + count: 1945560 + probability: 0.12975810913150784 +- from: 0 + to: 2 + count: 26381 + probability: 0.0017594670310852958 +- from: 0 + to: 3 + count: 2792 + probability: 0.000186210983313375 +- from: 0 + to: 4 + count: 0 + probability: 0.0 +- from: 1 + to: 0 + count: 1945560 + probability: 0.9666005228573729 +- from: 1 + to: 1 + count: 66078 + probability: 0.03282912341401421 +- from: 1 + to: 2 + count: 1123 + probability: 0.0005579331334776772 +- from: 1 + to: 3 + count: 25 + probability: 0.000012420595135300028 +- from: 1 + to: 4 + count: 0 + probability: 0.0 +- from: 2 + to: 0 + count: 26381 + probability: 0.9567346050627402 +- from: 2 + to: 1 + count: 1123 + probability: 0.04072677159643142 +- from: 2 + to: 2 + count: 51 + probability: 0.0018495684340320592 +- from: 2 + to: 3 + count: 19 + probability: 0.0006890549067962574 +- from: 2 + to: 4 + count: 0 + probability: 0.0 +- from: 3 + to: 0 + count: 2792 + probability: 0.9841381741275996 +- from: 3 + to: 1 + count: 25 + probability: 0.008812125484666901 +- from: 3 + to: 2 + count: 19 + probability: 0.0066972153683468455 +- from: 3 + to: 3 + count: 1 + probability: 0.00035248501938667606 +- from: 3 + to: 4 + count: 0 + probability: 0.0 +- from: 4 + to: 0 + count: 0 + probability: 0.0 +- from: 4 + to: 1 + count: 0 + probability: 0.0 +- from: 4 + to: 2 + count: 0 + probability: 0.0 +- from: 4 + to: 3 + count: 0 + probability: 0.0 +- from: 4 + to: 4 + count: 0 + probability: 0.0 +composition_transitions: +- from: 'A' + to: 'A' + count: 0 + probability: 0.0 +- from: 'A' + to: 'C' + count: 0 + probability: 0.0 +- from: 'A' + to: 'G' + count: 0 + probability: 0.0 +- from: 'A' + to: 'T' + count: 0 + probability: 0.0 +- from: 'C' + to: 'A' + count: 0 + probability: 0.0 +- from: 'C' + to: 'C' + count: 10771 + probability: 0.9619540948468339 +- from: 'C' + to: 'G' + count: 163 + probability: 0.014557470751094042 +- from: 'C' + to: 'T' + count: 263 + probability: 0.023488434402071982 +- from: 'G' + to: 'A' + count: 0 + probability: 0.0 +- from: 'G' + to: 'C' + count: 163 + probability: 0.006662851536952256 +- from: 'G' + to: 'G' + count: 24073 + probability: 0.9840173315892741 +- from: 'G' + to: 'T' + count: 228 + probability: 0.009319816873773708 +- from: 'T' + to: 'A' + count: 0 + probability: 0.0 +- from: 'T' + to: 'C' + count: 263 + probability: 0.008464484567603231 +- from: 'T' + to: 'G' + count: 228 + probability: 0.007338032248720672 +- from: 'T' + to: 'T' + count: 30580 + probability: 0.9841974831836761 diff --git a/doc/usage/phylo/index.html b/doc/usage/phylo/index.html index e19e2aca..c5488ce8 100644 --- a/doc/usage/phylo/index.html +++ b/doc/usage/phylo/index.html @@ -946,6 +946,17 @@ +
    • + +
  • + + + + Rare states and --iqtree-min-freq + + + +
    • @@ -1510,6 +1521,17 @@ + + +
  • + + + + Rare states and --iqtree-min-freq + + + +
    • @@ -1790,6 +1812,11 @@ Also write an IQ-TREE custom model and alignment (implies --sankoff) +--iqtree-min-freq +0.001 +With --iqtree --free-loss: also treat as missing any state rarer than this in the alignment + + --sankoff-cost-scale 100 Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands) @@ -2051,10 +2078,12 @@

    _sankoff.tcm (--phyg) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state 0, since it's never actually triggered).

    _sankoff.pg (--phyg) — script: read(prefasta:..., tcm:...) against _sankoff.fasta/_sankoff.tcm, a default 300s/4-instance search, report(...) writing <prefix>_sankoff.tre (bare numeric labels, as for --tnt). Run with phyg from the output directory (the script uses relative file names). Feed the tree plus _sankoff.fasta into name-tree for taxon names.

    _iqtree.model (--iqtree) — lower-triangular exchangeability matrix \(R(a,b) = e^{-\text{cost}(a,b)}\) (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered 0..k-1.

    -

    _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: +

    _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); with --iqtree-min-freq also set (the default), any state rarer than that threshold is folded into the same ? treatment, and non-informative columns are re-checked and dropped again after that. 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.

    +

    _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 — nor does any state --iqtree-min-freq folded away for being too rare. 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.

    +

    Rare states and --iqtree-min-freq

    +

    States that combine 3 or 4 central bases at once (IUPAC V/H/K.../N) are inherently rare — and, on real data, rare enough that they can make iqtree3 itself numerically unstable ("Numerical underflow for lh-derivative" warnings, near-degenerate likelihood optimization). They're also more likely to be assembly/detection noise than genuine, widely-shared multi-way polymorphism, the same "sampling failure, not true signal" reasoning --free-loss already applies to non-detection. With --free-loss set, --iqtree-min-freq (default 0.001, i.e. one in a thousand) extends that same missing-data treatment to any state below this frequency, not just absence. Check _iqtree_states.csv to see exactly which states survived and at what frequency; set --iqtree-min-freq 0 to disable this and keep every state that occurs at all (the old behavior). Has no effect without --free-loss — there is no missing-data symbol to fold rare states into otherwise.

    diff --git a/src/obikmer/src/cmd/phylo/args.rs b/src/obikmer/src/cmd/phylo/args.rs index 5261258e..32468a3c 100644 --- a/src/obikmer/src/cmd/phylo/args.rs +++ b/src/obikmer/src/cmd/phylo/args.rs @@ -267,6 +267,25 @@ pub struct PhyloArgs { #[arg(long)] pub iqtree: bool, + /// Under `--iqtree --free-loss`, also recode to `?` (the same + /// missing-data treatment as `-`) any state whose empirical frequency + /// in the alignment falls below this threshold — not just genuinely + /// absent calls. States encoding 3 or 4 simultaneously-observed + /// central bases (IUPAC `V`/`H`/`K`.../`N` for 3, `N` for 4) are rare + /// by construction and, on real data, land in exactly this low-frequency + /// range — more likely assembly/detection noise than a genuine, + /// widely-preserved multi-way polymorphism, the same "sampling failure, + /// not true signal" reasoning `--free-loss` already applies to absence. + /// Confirmed on real data to matter: `iqtree3`'s own "Numerical + /// underflow for lh-derivative" warnings and the exact-zero + /// exchangeability rows this was meant to fix (see + /// `DevDocMD/theory/evolutionary_distances.md`) both trace back to + /// states this thin. No effect without `--free-loss` (there is no + /// missing-data symbol to recode to otherwise). `_iqtree_states.csv` + /// reports the frequency actually used to decide. + #[arg(long, default_value = "0.001")] + pub iqtree_min_freq: f64, + /// Scale factor applied before rounding real-valued costs to the /// integers both `--tnt`'s smatrix/cost commands and `--phyg`'s `tcm:` /// matrix require. Keep this small: the total tree score is this scale diff --git a/src/obikmer/src/cmd/phylo/iqtree.rs b/src/obikmer/src/cmd/phylo/iqtree.rs index b774ce59..3ec1b712 100644 --- a/src/obikmer/src/cmd/phylo/iqtree.rs +++ b/src/obikmer/src/cmd/phylo/iqtree.rs @@ -107,6 +107,19 @@ fn drop_ascertainment_noninformative(alignment: &SnpAlignment) -> SnpAlignment { SnpAlignment { sequences } } +/// Recode every occurrence of a byte in `symbols` to `-` — the same +/// "absent" byte `drop_ascertainment_noninformative`/`compact_alphabet` +/// already treat specially under `--free-loss` (recoded to `?` further +/// downstream). Used by `--iqtree-min-freq` to fold rare, likely-noisy +/// states into the missing-data treatment before a second +/// `compact_alphabet` pass, without duplicating that treatment's logic. +fn recode_symbols_as_absent(alignment: &SnpAlignment, symbols: &[u8]) -> SnpAlignment { + let sequences = alignment.sequences.iter() + .map(|seq| seq.iter().map(|&b| if symbols.contains(&b) { b'-' } else { b }).collect()) + .collect(); + SnpAlignment { sequences } +} + fn compact_alphabet(alignment: &SnpAlignment, free_loss: bool) -> CompactAlphabet { let iupac_to_state = state_index_table(); @@ -240,6 +253,7 @@ pub(super) fn write_iqtree( labels: &[String], output: &Option, free_loss: bool, + min_freq: f64, ) { let filtered; let alignment = if free_loss { @@ -257,7 +271,41 @@ pub(super) fn write_iqtree( alignment }; - let alphabet = compact_alphabet(alignment, free_loss); + let mut alphabet = compact_alphabet(alignment, free_loss); + + // `--iqtree-min-freq`: fold rare (likely-noisy) states into the same + // missing-data treatment `-` already gets under `--free-loss`, then + // recompute the alphabet on the further-filtered alignment — see + // `args.rs`'s docs on `--iqtree-min-freq` and + // `DevDocMD/theory/evolutionary_distances.md` for why (real data: + // states encoding 3-4 simultaneous central bases land in exactly this + // low-frequency range and correlate with `iqtree3`'s own numerical + // instability warnings). + let refiltered; + let alignment = if free_loss { + let low_freq_symbols: Vec = alphabet.compact_to_old.iter().zip(alphabet.freq.iter()) + .filter(|&(_, &f)| f < min_freq) + .map(|(&old, _)| STATE_SYMBOL[old as usize] as u8) + .collect(); + if low_freq_symbols.is_empty() { + alignment + } else { + let recoded = recode_symbols_as_absent(alignment, &low_freq_symbols); + let before = recoded.sequences.first().map(|s| s.len()).unwrap_or(0); + refiltered = drop_ascertainment_noninformative(&recoded); + let after = refiltered.sequences.first().map(|s| s.len()).unwrap_or(0); + info!( + "--iqtree-min-freq {min_freq}: {} rare state(s) ({}) recoded as missing, {before} → {after} sites", + low_freq_symbols.len(), + low_freq_symbols.iter().map(|&b| b as char).collect::(), + ); + alphabet = compact_alphabet(&refiltered, free_loss); + &refiltered + } + } else { + alignment + }; + 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); @@ -349,4 +397,39 @@ mod tests { assert_eq!(lines.next(), Some("1,G,0.5")); assert!(lines.next().is_none()); } + + #[test] + fn iqtree_min_freq_folds_rare_states_into_missing() { + // 20 common A/C sites (60 calls total across 3 genomes) plus one + // site where genome 0 carries the rare ambiguity state `M` and + // genome 1 carries `A` (kept informative by the first + // ascertainment filter: two distinct non-`-` calls) — `M` ends up + // at 1/62, well below the 0.05 threshold used here. + let mut sequences: Vec> = vec![Vec::new(); 3]; + for i in 0..20 { + let (a, b, c) = if i % 2 == 0 { (b'A', b'C', b'A') } else { (b'C', b'A', b'C') }; + sequences[0].push(a); + sequences[1].push(b); + sequences[2].push(c); + } + sequences[0].push(b'M'); + sequences[1].push(b'A'); + sequences[2].push(b'-'); + let alignment = SnpAlignment { sequences }; + let labels = vec!["g1".to_string(), "g2".to_string(), "g3".to_string()]; + let matrix = [[0.0f64; 16]; 16]; + let prefix = std::env::temp_dir().join(format!("obikmer_test_iqtree_minfreq_{}", std::process::id())); + let output = Some(prefix.clone()); + + write_iqtree(&matrix, &alignment, &labels, &output, true, 0.05); + + let states_path = format!("{}_iqtree_states.csv", prefix.display()); + let csv = std::fs::read_to_string(&states_path).unwrap(); + assert!(!csv.contains(",M,"), "M (freq ~1/62) must be folded into missing under --iqtree-min-freq 0.05, got:\n{csv}"); + assert!(csv.contains(",A,") && csv.contains(",C,"), "A/C must survive (well above threshold), got:\n{csv}"); + + for suffix in ["_iqtree_states.csv", "_iqtree.model", "_iqtree.fasta"] { + std::fs::remove_file(format!("{}{suffix}", prefix.display())).ok(); + } + } } diff --git a/src/obikmer/src/cmd/phylo/mod.rs b/src/obikmer/src/cmd/phylo/mod.rs index b3684ccf..cfeaf816 100644 --- a/src/obikmer/src/cmd/phylo/mod.rs +++ b/src/obikmer/src/cmd/phylo/mod.rs @@ -249,7 +249,7 @@ pub fn run(args: PhyloArgs) { write_sankoff_phyg(&matrix, &args.output, args.sankoff_cost_scale); } if args.iqtree { - write_iqtree(&matrix, &alignment, &kept_labels, &args.output, args.free_loss); + write_iqtree(&matrix, &alignment, &kept_labels, &args.output, args.free_loss, args.iqtree_min_freq); } } diff --git a/src/obikphylo/src/siblings/tests.rs b/src/obikphylo/src/siblings/tests.rs index 2f30edaa..f984a376 100644 --- a/src/obikphylo/src/siblings/tests.rs +++ b/src/obikphylo/src/siblings/tests.rs @@ -15,7 +15,7 @@ use super::cardinality::CardinalityExt; use super::distance::DistanceExt; use super::entropy::ShannonEntropyExt; use super::entropy_annex::{EntropyAnnex, ENTROPY_ANNEX_FILE_NAME}; -use super::helpers::is_minorant; +use super::helpers::{central_base, is_minorant}; use super::sankoff_bundle::SankoffBundleExt; use super::stats::SiblingStatsExt; use super::subsample::EntropyBias; @@ -457,6 +457,44 @@ fn sankoff_bundle_matches_old_separate_calls() { assert_eq!(bundle.alignment.sequences, expected_alignment.sequences); } +#[test] +fn base_pair_tally_accumulates_base_a_diagnostic() { + // Diagnostic for a user-reported observation on real data: + // `composition_transitions`'s row/column for base 'A' (index 0) was + // *entirely* zero, including the diagonal (`same[0]`, "A stays A"), + // despite 'A' appearing at ~6% frequency in the pseudo-alignment + // itself. g1/g2 are identical (both central base 'A'), g3 differs + // (central base 'C') at the same family — makes the family variable + // (family_size=2) while g1/g2 form an eligible (ratio 0), ordinary + // pair both resolving `single_form` to `Some(0)` = 'A'. + let dir = tempdir().unwrap(); + let g1 = build_single_genome_index(dir.path(), "g1", b"AAAAAAAAAAA"); + let g2 = build_single_genome_index(dir.path(), "g2", b"AAAAAAAAAAA"); + let g3 = build_single_genome_index(dir.path(), "g3", b"AAAAACAAAAA"); + let mut rep = Reporter::new(); + let merged = KmerIndex::merge( + &dir.path().join("merged.idx"), + &[&g1, &g2, &g3], + MergeMode::Presence, + false, + false, + 1.0, + &mut rep, + ).expect("merge"); + merged.build_sibling_annex().expect("build_sibling_annex"); + + let n_genomes = merged.meta().genomes.len(); + let exclude_mask = vec![false; n_genomes]; + let bundle = merged.sankoff_bundle(None, None, 0.5, &exclude_mask).expect("sankoff_bundle"); + + assert!( + bundle.base_pair_tally.same[0] >= 1, + "base 'A' (index 0) must accumulate in BasePairTally.same for an identical, included g1/g2 \ + pair both resolving to 'A' — got same={:?}, counts={:?}", + bundle.base_pair_tally.same, bundle.base_pair_tally.counts, + ); +} + #[test] fn entropy_annex_builds_on_demand_and_biases_selection() { // Same fixture, same single non-monomorphic family, known entropy15 = @@ -537,6 +575,132 @@ fn diag_real_index_layer_distribution() { println!("top10 counts: {:?}", &sorted_counts[..10.min(sorted_counts.len())]); } +#[test] +#[ignore] +fn diag_real_index_base_a_representation() { + // Diagnostic for a user-reported observation, reproduced on two + // unrelated real datasets (a plant and this bacterial index): + // `composition_transitions`'s row/column for base 'A' is entirely + // zero, even though 'A' appears at real (~6-22%) frequency in the + // pseudo-alignment itself. A minimal synthetic reproduction + // (`base_pair_tally_accumulates_base_a_diagnostic`) showed the + // pairwise tally mechanism itself works correctly for base A in + // isolation — so this checks one level lower, purely structural + // (annex bits only, no cross-partition/genome-level resolution): does + // base A even show up as *present* in minorant families' own + // `FamilyMask`s at all, at a rate proportional to the other 3 bases? + // If yes, the anomaly is specific to the pairwise/`single_form` + // resolution stage; if `has_base[0]` is itself near-zero here, the + // anomaly originates earlier, in `build_sibling_annex`/`central_base`. + let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence") + .expect("open real index"); + let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs"); + let mut has_base = [0u64; 4]; + let mut minorant_total = 0u64; + let mut minorant_central_a = 0u64; + for layer_dir in &layer_dirs { + let annex = SiblingAnnex::open(&layer_dir.join(ANNEX_FILE_NAME)).unwrap(); + for slot in 0..annex.len() { + let Some(mask) = annex.get(slot) else { continue }; + if !mask.is_minorant() { + continue; + } + minorant_total += 1; + for b in 0..4u8 { + if mask.has(b) { + has_base[b as usize] += 1; + } + } + } + // Cross-check: for a sample of this layer's minorants, is the + // slot's *own* central base ever actually 'A' (bit 0)? — reads + // the real k-mer via the layer's MPHF, not just the mask. + let meta = PartitionMeta::load(layer_dir.parent().unwrap()).unwrap(); + let mphf = MphfLayer::open(layer_dir, &meta.mode).unwrap(); + for (order, kmer) in mphf.enumerate_kmers().take(2_000_000) { + let Some(mask) = annex.get(order) else { continue }; + if mask.is_minorant() && central_base(kmer, idx.kmer_size()) == 0 { + minorant_central_a += 1; + } + } + } + println!( + "minorant_total={minorant_total} has_base(A,C,G,T)={has_base:?} minorant_central_a_sampled={minorant_central_a}" + ); +} + +#[test] +#[ignore] +fn diag_real_index_sankoff_bundle_base_a() { + // Structural check (`diag_real_index_base_a_representation`) shows + // base A fully, proportionally represented at the family level (~23M + // minorants, the *highest* of the four bases) — so the anomaly must + // be downstream, in `sankoff_bundle`'s actual pairwise resolution. + // Reproduces through the real cross-partition code path (not the + // minimal synthetic fixture, which showed the mechanism working in + // isolation), bounded by `--subsample` to stay fast. + let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence") + .expect("open real index"); + let n_genomes = idx.meta().genomes.len(); + let exclude_mask = vec![false; n_genomes]; + let bundle = idx.sankoff_bundle(Some(1_000_000), None, 0.5, &exclude_mask).expect("sankoff_bundle"); + println!("base_pair_tally.same={:?}", bundle.base_pair_tally.same); + println!("base_pair_tally.counts={:?}", bundle.base_pair_tally.counts); + let alignment_a_count: usize = bundle.alignment.sequences.iter() + .map(|seq| seq.iter().filter(|&&b| b == b'A').count()) + .sum(); + println!("alignment raw 'A' byte count={alignment_a_count}"); +} + +#[test] +#[ignore] +fn diag_real_index_genome_mask_for_base_a_families() { + // Directly inspects the raw `genome_mask` array `sankoff_bundle`'s + // closures see, for the first few variable families where base A is + // present — to check whether A ever co-occurs as `single_form` in two + // *different* genomes at the same family at all (the precondition for + // `bp_same`/`bp_counts` to ever increment for base A), rather than + // reasoning about it further. + use std::sync::Arc; + let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence") + .expect("open real index"); + let n_parts = idx.n_partitions(); + let n_genomes = idx.meta().genomes.len(); + let with_counts = idx.meta().config.with_counts; + let k = idx.kmer_size(); + let n_bits = n_parts.trailing_zeros() as usize; + let partition = obikpartitionner::KmerPartition::open_with_config( + idx.root_path(), idx.kmer_size(), idx.minimizer_size(), n_bits, + ).unwrap(); + let cache = Arc::new(super::cache::PartitionCache::build(&partition, n_parts, with_counts).unwrap()); + let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).unwrap(); + + let mut printed = 0usize; + for layer_dir in &layer_dirs { + super::family_scan::scan_layer_families( + layer_dir, n_parts, n_genomes, with_counts, k, &cache, &super::family_scan::Selection::All, + |_family_idx, mask, genome_mask| { + if printed >= 15 || !mask.has(0) || mask.family_size() < 2 { + return; + } + let single_a: Vec = (0..n_genomes).filter(|&g| genome_mask[g] == 1).collect(); + let others: Vec<(usize, u8)> = (0..n_genomes) + .filter(|&g| genome_mask[g] != 0 && genome_mask[g] != 1) + .map(|g| (g, genome_mask[g])) + .collect(); + println!( + "family bits={:#06b} size={} single_A_genomes={:?} other_nonzero_genomes={:?}", + mask.bits(), mask.family_size(), single_a, others, + ); + printed += 1; + }, + ).unwrap(); + if printed >= 15 { + break; + } + } +} + #[test] #[ignore] fn diag_plant_index_presence_matrix_sparsity() {