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.
This commit is contained in:
@@ -2051,13 +2051,28 @@ but indexed by <code>family_idx</code> (every minorant of the layer, monomorphic
|
|||||||
included — the same numbering <code>Selection</code>/<code>scan_layer_families</code> already
|
included — the same numbering <code>Selection</code>/<code>scan_layer_families</code> already
|
||||||
use), one <code>f32</code> entropy15 value per entry, <code>-1.0</code> sentinel for monomorphic/
|
use), one <code>f32</code> entropy15 value per entry, <code>-1.0</code> sentinel for monomorphic/
|
||||||
not-yet-computed. First use of <code>--entropy</code>/<code>--entropy-sd</code> on an index
|
not-yet-computed. First use of <code>--entropy</code>/<code>--entropy-sd</code> on an index
|
||||||
pays a one-time cost (<code>ensure_entropy_annexes</code> in <code>entropy.rs</code>: a full,
|
pays a one-time cost (<code>ensure_entropy_annexes</code> in <code>entropy.rs</code>) — every
|
||||||
unsampled <code>Selection::All</code> scan, resolving every non-monomorphic
|
|
||||||
minorant's <code>genome_mask</code> once to compute and persist its entropy) — every
|
|
||||||
later run (any <code>μ</code>/<code>σ</code>, any command) reads the file positionally, no
|
later run (any <code>μ</code>/<code>σ</code>, any command) reads the file positionally, no
|
||||||
re-scan, restoring the usual <code>Selection::Some</code> "skip resolving excluded
|
re-scan, restoring the usual <code>Selection::Some</code> "skip resolving excluded
|
||||||
families" speedup that a naive "weigh during the resolving scan" design
|
families" speedup that a naive "weigh during the resolving scan" design
|
||||||
would have permanently forfeited.</p>
|
would have permanently forfeited.</p>
|
||||||
|
<p><strong>Bug found and fixed (2026-08-15): <code>ensure_entropy_annexes</code> scanned with
|
||||||
|
<code>Selection::All</code> instead of bounding to non-monomorphic minorants.</strong>
|
||||||
|
Monomorphism (<code>family_size() < 2</code>) is knowable directly from the annex
|
||||||
|
bits alone, no per-genome resolution needed — but the original
|
||||||
|
implementation called <code>scan_layer_families</code> with <code>Selection::All</code>
|
||||||
|
anyway, so <code>fill_sub_matrix_carries</code> (the expensive per-genome
|
||||||
|
resolution) ran for <em>every</em> minorant, ~98% of which are monomorphic
|
||||||
|
(measured elsewhere in this doc) and had their <code>genome_mask</code> immediately
|
||||||
|
discarded once the callback checked <code>family_size() < 2</code>. Fixed by adding
|
||||||
|
[<code>subsample::non_monomorphic_selection_layer</code>] — a cheap, annex-only,
|
||||||
|
non-sampling pass (same shape as <code>reservoir_sample_layer</code>, but keeping
|
||||||
|
every non-monomorphic minorant's <code>family_idx</code> instead of a bounded
|
||||||
|
reservoir) — and passing <code>Selection::Some(&eligible)</code> instead of
|
||||||
|
<code>Selection::All</code>, so the expensive resolution now runs only for the ~2%
|
||||||
|
of minorants that can actually produce a real entropy value. A
|
||||||
|
<code>debug_assert!(mask.family_size() >= 2, ...)</code> inside the
|
||||||
|
<code>scan_layer_families</code> callback guards the invariant.</p>
|
||||||
<p><strong>Resolved</strong>: the existing hard "non-monomorphic minorant" eligibility
|
<p><strong>Resolved</strong>: the existing hard "non-monomorphic minorant" eligibility
|
||||||
filter stays a hard gate upstream of the Gaussian weighting — only
|
filter stays a hard gate upstream of the Gaussian weighting — only
|
||||||
qualifying families ever get a stored entropy value or a weighted draw.</p>
|
qualifying families ever get a stored entropy value or a weighted draw.</p>
|
||||||
|
|||||||
@@ -887,6 +887,17 @@
|
|||||||
</ul>
|
</ul>
|
||||||
</nav>
|
</nav>
|
||||||
|
|
||||||
|
</li>
|
||||||
|
|
||||||
|
<li class="md-nav__item">
|
||||||
|
<a href="#_iqtree_statescsv-compact-symbol-traceability-2026-08-15" class="md-nav__link">
|
||||||
|
<span class="md-ellipsis">
|
||||||
|
|
||||||
|
_iqtree_states.csv: compact-symbol traceability (2026-08-15)
|
||||||
|
|
||||||
|
</span>
|
||||||
|
</a>
|
||||||
|
|
||||||
</li>
|
</li>
|
||||||
|
|
||||||
<li class="md-nav__item">
|
<li class="md-nav__item">
|
||||||
@@ -2095,6 +2106,17 @@
|
|||||||
</ul>
|
</ul>
|
||||||
</nav>
|
</nav>
|
||||||
|
|
||||||
|
</li>
|
||||||
|
|
||||||
|
<li class="md-nav__item">
|
||||||
|
<a href="#_iqtree_statescsv-compact-symbol-traceability-2026-08-15" class="md-nav__link">
|
||||||
|
<span class="md-ellipsis">
|
||||||
|
|
||||||
|
_iqtree_states.csv: compact-symbol traceability (2026-08-15)
|
||||||
|
|
||||||
|
</span>
|
||||||
|
</a>
|
||||||
|
|
||||||
</li>
|
</li>
|
||||||
|
|
||||||
<li class="md-nav__item">
|
<li class="md-nav__item">
|
||||||
@@ -4048,6 +4070,42 @@ with a different taxon count crashes with <code>Assertion 'i >= 0 &&
|
|||||||
the old leaf count onto the new, smaller alignment. Not an <code>obikmer</code> bug;
|
the old leaf count onto the new, smaller alignment. Not an <code>obikmer</code> bug;
|
||||||
avoid by using a fresh <code>--prefix</code> (or <code>-redo</code>) whenever the underlying
|
avoid by using a fresh <code>--prefix</code> (or <code>-redo</code>) whenever the underlying
|
||||||
alignment's taxon set changes, never <code>--undo</code> across them.</p>
|
alignment's taxon set changes, never <code>--undo</code> across them.</p>
|
||||||
|
<h2 id="_iqtree_statescsv-compact-symbol-traceability-2026-08-15"><code>_iqtree_states.csv</code>: compact-symbol traceability (2026-08-15)</h2>
|
||||||
|
<p>User-reported suspicion: <code>--iqtree --free-loss</code> state frequencies "don't
|
||||||
|
sum to 1". Verified both by code trace and unit test
|
||||||
|
(<code>obikmer/src/cmd/phylo/iqtree.rs</code>, <code>free_loss_excludes_absent_state_and_freq_sums_to_one</code>)
|
||||||
|
— <strong>not a bug</strong>: <code>compact_alphabet</code>'s counting loop <code>continue</code>s on
|
||||||
|
<code>free_loss && b == b'-'</code> <em>before</em> recoding/counting, so the absent state
|
||||||
|
never enters <code>occurs</code>/<code>counts</code>, and <code>freq[i] = counts[old]/total</code> sums to
|
||||||
|
1.0 by construction (<code>total</code> is defined as the sum over exactly the same
|
||||||
|
states). Confirmed against real data the user provided: a pasted
|
||||||
|
<code>_iqtree.model</code> frequency line summed to <code>1.000001</code> — a 6-decimal
|
||||||
|
print-rounding artifact (<code>format!("{p:.6}")</code> × 15 values), not a
|
||||||
|
computation error; IQ-TREE's own 3-decimal-rounded report of the same
|
||||||
|
values summed to exactly <code>1.000</code>.</p>
|
||||||
|
<p>Investigating that report surfaced a real, separate gap: nothing mapped
|
||||||
|
<code>_iqtree.model</code>/<code>_iqtree.fasta</code>'s compact state symbols (<code>0-9A-F</code>, IQ-TREE
|
||||||
|
renumbers away unused states from the full 16, see this module's own
|
||||||
|
"<code>--iqtree</code>... state count" discussion) back to the canonical 16-state
|
||||||
|
alphabet (<code>STATE_SYMBOL</code>, <code>sankoff.rs</code>) — so a pattern like "compact state
|
||||||
|
0 has zero exchangeability with every other state" (<code>R(a,b) =
|
||||||
|
exp(-cost(a,b)) = 0</code> 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 <code>cost = -ln(0) = +∞</code>) or a genuine problem.</p>
|
||||||
|
<p><strong>Fix</strong>: <code>write_iqtree_states_csv</code> (<code>iqtree.rs</code>) writes
|
||||||
|
<code><prefix>_iqtree_states.csv</code> — one row per surviving state,
|
||||||
|
<code>iqtree_symbol,canonical_symbol,frequency</code>, frequency at full <code>f64</code>
|
||||||
|
precision (not the model file's truncated 6 decimals). Written alongside
|
||||||
|
<code>_iqtree.model</code>/<code>_iqtree.fasta</code> from the same <code>CompactAlphabet</code> both
|
||||||
|
already use, so there is no risk of the three files disagreeing. Covered
|
||||||
|
by <code>states_csv_maps_compact_symbols_back_to_canonical_ones</code>.</p>
|
||||||
|
<p>The zero-exchangeability pattern itself (state 0 in the user's report,
|
||||||
|
frequency 3.26%, <code>R=0</code> with every other state) is not yet explained —
|
||||||
|
plausibly a genuinely unobserved transition in the calibration
|
||||||
|
(<code>cardinality_transitions</code>/<code>composition_transitions</code> count <code>0</code> for every
|
||||||
|
pair involving it), which is a legitimate, if numerically extreme, result
|
||||||
|
of <code>-ln(0)</code>, not necessarily a bug — not investigated further.</p>
|
||||||
<h2 id="references">References</h2>
|
<h2 id="references">References</h2>
|
||||||
<p>The Mash mutation-rate model this discussion contrasts with:
|
<p>The Mash mutation-rate model this discussion contrasts with:
|
||||||
(Fan <em>et al.</em> 2015; Marbl Lab 2026)<sup id="fnref:Mash-distances-doc"><a class="footnote-ref" href="#fn:Mash-distances-doc">1</a></sup> <sup id="fnref:Fan2015-mash-formula"><a class="footnote-ref" href="#fn:Fan2015-mash-formula">2</a></sup>.</p>
|
(Fan <em>et al.</em> 2015; Marbl Lab 2026)<sup id="fnref:Mash-distances-doc"><a class="footnote-ref" href="#fn:Mash-distances-doc">1</a></sup> <sup id="fnref:Fan2015-mash-formula"><a class="footnote-ref" href="#fn:Fan2015-mash-formula">2</a></sup>.</p>
|
||||||
|
|||||||
@@ -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
|
avoid by using a fresh `--prefix` (or `-redo`) whenever the underlying
|
||||||
alignment's taxon set changes, never `--undo` across them.
|
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
|
||||||
|
`<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
|
## References
|
||||||
|
|
||||||
The Mash mutation-rate model this discussion contrasts with:
|
The Mash mutation-rate model this discussion contrasts with:
|
||||||
|
|||||||
@@ -210,6 +210,7 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below):
|
|||||||
| `<prefix>_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search |
|
| `<prefix>_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search |
|
||||||
| `<prefix>_iqtree.model` | `--iqtree` | IQ-TREE model file | custom ML substitution model |
|
| `<prefix>_iqtree.model` | `--iqtree` | IQ-TREE model file | custom ML substitution model |
|
||||||
| `<prefix>_iqtree.fasta` | `--iqtree` | FASTA | alignment recoded for that model |
|
| `<prefix>_iqtree.fasta` | `--iqtree` | 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.
|
**`_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 <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
|
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.
|
||||||
|
|||||||
@@ -2014,6 +2014,12 @@
|
|||||||
<td>FASTA</td>
|
<td>FASTA</td>
|
||||||
<td>alignment recoded for that model</td>
|
<td>alignment recoded for that model</td>
|
||||||
</tr>
|
</tr>
|
||||||
|
<tr>
|
||||||
|
<td><code><prefix>_iqtree_states.csv</code></td>
|
||||||
|
<td><code>--iqtree</code></td>
|
||||||
|
<td>CSV table</td>
|
||||||
|
<td>maps <code>_iqtree.model</code>/<code>_iqtree.fasta</code>'s compact state symbols back to <code>_sankoff_matrix.csv</code>'s alphabet</td>
|
||||||
|
</tr>
|
||||||
</tbody>
|
</tbody>
|
||||||
</table>
|
</table>
|
||||||
<p><strong><code>_sankoff_matrix.csv</code></strong> — header <code>state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N</code>: 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), <code>0</code> standing for the empty/absent state (not <code>-</code>, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost <span class="arithmatex">\(-\ln P(a,b)\)</span>, 4 decimals.</p>
|
<p><strong><code>_sankoff_matrix.csv</code></strong> — header <code>state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N</code>: 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), <code>0</code> standing for the empty/absent state (not <code>-</code>, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost <span class="arithmatex">\(-\ln P(a,b)\)</span>, 4 decimals.</p>
|
||||||
@@ -2048,6 +2054,7 @@
|
|||||||
<p><strong><code>_iqtree.fasta</code></strong> (<code>--iqtree</code>) — alignment recoded to that same compact <code>0..k-1</code> alphabet (symbols <code>0-9A-F</code>). Under <code>--free-loss</code>, non-detection becomes <code>?</code> and columns left non-informative once missing calls are ignored are dropped first (required for <code>+ASC</code>). Run with:
|
<p><strong><code>_iqtree.fasta</code></strong> (<code>--iqtree</code>) — alignment recoded to that same compact <code>0..k-1</code> alphabet (symbols <code>0-9A-F</code>). Under <code>--free-loss</code>, non-detection becomes <code>?</code> and columns left non-informative once missing calls are ignored are dropped first (required for <code>+ASC</code>). Run with:
|
||||||
<div class="highlight"><pre><span></span><code>iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
|
<div class="highlight"><pre><span></span><code>iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
|
||||||
</code></pre></div></p>
|
</code></pre></div></p>
|
||||||
|
<p><strong><code>_iqtree_states.csv</code></strong> (<code>--iqtree</code>) — one row per state actually kept in <code>_iqtree.model</code>/<code>_iqtree.fasta</code> (header <code>iqtree_symbol,canonical_symbol,frequency</code>): <code>iqtree_symbol</code> is the compact <code>0-9A-F</code> symbol as written in those two files, <code>canonical_symbol</code> is the matching <code>_sankoff_matrix.csv</code> state, <code>frequency</code> is that state's empirical frequency at full precision (<code>_iqtree.model</code>'s own frequency line is rounded to 6 decimals). Under <code>--free-loss</code>, absent (<code>0</code>/<code>?</code>) is never a kept state, so it never appears here. Use this file to identify which real state a given row/column of <code>_iqtree.model</code>'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.</p>
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@@ -5,7 +5,7 @@ use obifastwrite::{JsonVal, write_record};
|
|||||||
use obikphylo::siblings::SnpAlignment;
|
use obikphylo::siblings::SnpAlignment;
|
||||||
use tracing::info;
|
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 ──
|
// ── 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 }
|
CompactAlphabet { old_to_compact, compact_to_old, freq }
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/// Write `<prefix>_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<PathBuf>) -> 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
|
/// Write the `R` (exchangeability) + `π` (frequencies) model file IQ-TREE's
|
||||||
/// `-m <file>+ASC` reads. Returns the path, so the caller can print a
|
/// `-m <file>+ASC` reads. Returns the path, so the caller can print a
|
||||||
/// single combined "how to run this" message once the alignment is also
|
/// 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 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 model_path = write_iqtree_model(matrix, &alphabet, output);
|
||||||
let (fasta_path, n_sites) = write_iqtree_alignment(alignment, labels, &alphabet, output, free_loss);
|
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());
|
.unwrap_or_else(|| "iqtree".into());
|
||||||
info!(
|
info!(
|
||||||
"IQ-TREE alignment → {fasta_path} ({n_sites} sites, {} states)\n\
|
"IQ-TREE alignment → {fasta_path} ({n_sites} sites, {} states)\n\
|
||||||
|
IQ-TREE state mapping → {states_path}\n\
|
||||||
Run with:\n \
|
Run with:\n \
|
||||||
iqtree3 -s {fasta_path} --seqtype MORPH -m {model_path}+ASC --prefix {prefix_name} -T AUTO",
|
iqtree3 -s {fasta_path} --seqtype MORPH -m {model_path}+ASC --prefix {prefix_name} -T AUTO",
|
||||||
alphabet.k()
|
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());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user