feat(distance): implement native Sankoff calibration and backends
Replaces external Python glue with native Rust modules for Sankoff model calibration, exporting calibrated cost matrices, FASTA alignments, and YAML parameters. Adds dedicated writers for TNT and PhyG that apply integer scaling and Floyd-Warshall metric closure to enforce triangle inequality. Integrates these exporters into the distance command pipeline to streamline downstream tree inference workflows, while updating theory documentation to reflect IQ-TREE integration and state renumbering improvements.
This commit is contained in:
Eric Coissac
@@ -620,6 +620,285 @@ simply does not survive domain-level divergence. Cross-domain placement
|
||||
would need conserved-marker characters (rRNA, ribosomal proteins), not this
|
||||
estimator.
|
||||
|
||||
### Native `--sankoff --tnt`/`--phyg` export (2026-08-12), superseding the external scripts above
|
||||
|
||||
The ad hoc Python glue from the previous section is superseded: `obikmer
|
||||
distance --sankoff` now calibrates the matrix natively (`p_hat`, 6-category
|
||||
substitution costs, `c_ctx` weighted by `mean_sub_cost` — see the worked
|
||||
example above) and `--tnt`/`--phyg` each write a ready-to-run script from
|
||||
it, no external script needed. `∅` is an ordinary 16th state throughout
|
||||
(never `-`), specifically to avoid gap-semantics confusion in downstream
|
||||
tools — see "A concrete Sankoff cost matrix" above for why.
|
||||
|
||||
**TNT (`--tnt`).** `write_sankoff_tnt` (`obikmer/src/cmd/distance.rs`)
|
||||
recodes to TNT's own `0-9A-F` xread alphabet (its default reader rejects
|
||||
the wider IUPAC set otherwise), scales and rounds costs to integers
|
||||
(`smatrix`/`cost` reject decimals), then re-runs integer Floyd-Warshall on
|
||||
the rounded matrix (`scaled_metric_matrix`) — independently rounding each
|
||||
cell of an already-metric real-valued matrix can break the triangle
|
||||
inequality (e.g. two real costs of `1.734` round to `173` each, summing to
|
||||
`346`, while their real sum `3.468` rounds to `347`), which TNT otherwise
|
||||
silently "fixes" itself with an unreproducible correction. Verified against
|
||||
the real 20-genome benchmark index: zero triangle-inequality violations
|
||||
after the fix, TNT loads the file without its "triangle inequality
|
||||
violated... Fixed" warning.
|
||||
|
||||
Two syntax facts worth recording because they're wrong in intuitive
|
||||
guesses and contradicted actual TNT behavior when tested: TNT's plain
|
||||
command stream has **no comment syntax** of its own — `/* */` and `[ ]`
|
||||
only work inside the (separately-enabled, off by default) macro scripting
|
||||
language, and error with "No command!" otherwise. The working substitute
|
||||
is `quote TEXT ;` (prints the text, doesn't affect parsing) — but the text
|
||||
itself can't contain a literal `;` (TNT's universal terminator); the
|
||||
manual's own escape (`.,`) exists but the script here just avoids
|
||||
semicolons in the text instead.
|
||||
|
||||
The default search command embedded in the script is `mult` (traditional:
|
||||
random addition sequences + TBR), not `xmult` (New Technology search:
|
||||
ratchet/drift/tree-fusion). `xmult` with TNT's default `mxram` (16 MB,
|
||||
must be set *before* `xread` if changed) ran out of RAM on the real
|
||||
908k-character dataset ("`xmult - out of ram`"); `mult` does not, matching
|
||||
what had already been validated by hand outside this session.
|
||||
|
||||
**PhyG (`--phyg`).** `write_sankoff_phyg` writes a `tcm:` custom-alphabet
|
||||
matrix (same scale+round+metric-closure treatment as TNT) and reuses
|
||||
`--sankoff`'s own `_sankoff.fasta` as-is via `prefasta:` — PhyG's `tcm:`
|
||||
alphabet is read from the matrix file's own first line, so (unlike TNT) no
|
||||
recoding is needed. PhyG auto-adds its own indel/gap state as an
|
||||
`(n+1)`-th row/column of the tcm; inert here since the alignment encodes
|
||||
absence as `0`, never `-`.
|
||||
|
||||
`report("file", newick, overwrite)` — exactly as shown in PhyG's own
|
||||
manual — triggers `Unrecognized/missing report option ... defaulting to
|
||||
'graphs'` on the locally installed binary (1.3, commit `3c1a1fa`); the
|
||||
working form adds `graphs` explicitly: `report("file", graphs, newick,
|
||||
overwrite)`. Manual/binary mismatches like this (also true of `criterion:`
|
||||
— the binary accepts `parsimony`/`ml`/`pmdl`, the manual instead documents
|
||||
`mapa`/`ncm`/`parsimony`/`pmdl`/`si`) mean command syntax against this PhyG
|
||||
build should be verified empirically, not trusted from the PDF alone.
|
||||
`instances:N` (not a separate CPU flag) is what actually parallelises the
|
||||
search across cores — PhyG uses all physical cores by default but only
|
||||
across as many instances as are running, so raise it to the physical core
|
||||
count to use them all (the CLI-level `+RTS -NX -RTS` flag also exists but
|
||||
controls something else: capping/limiting cores, not requesting more).
|
||||
|
||||
Both scripts share one `--sankoff-cost-scale` (default `100`), not two
|
||||
separate flags — they scale the same calibrated matrix for the same
|
||||
reason (integer-only cost commands) and no PhyG-specific
|
||||
accumulator-width constraint was ever found to justify a different
|
||||
default from TNT's (TNT: hinted 32-bit accumulators in its own manual;
|
||||
PhyG: no such hint found — Haskell's native `Int` is typically 64-bit).
|
||||
|
||||
**Open problem: PhyG reports all branch lengths as `0.0`.** The graph-level
|
||||
parsimony cost is correct (`3.3286×10⁸` on the real dataset, consistent
|
||||
with TNT's `328574911` on the same calibrated matrix), but every individual
|
||||
edge in the exported Newick shows `:0.0`, with the total cost only ever
|
||||
shown as a whole-tree annotation (`[3.32860377e8]`). Not fixed, not fully
|
||||
diagnosed — PhyG's manual describes per-edge branch length as computed by
|
||||
ancestral-state (HTU) backtracking, well documented for sequence/standard
|
||||
character types, but nothing found (the term "Sankoff" doesn't even appear
|
||||
in the manual) confirming this backtracking is wired up for a custom
|
||||
`tcm:` matrix character. Switching `criterion:` to a likelihood-family
|
||||
option (`ml`, or the manual's `mapa`/`ncm`/`si`) was considered as a
|
||||
possible fix but is very unlikely to be one: those criteria are
|
||||
information-theoretic reparametrisations of the *same* step-counting
|
||||
machinery as parsimony (`ncm` in particular is known in the literature to
|
||||
be numerically equivalent to weighted parsimony), not classical
|
||||
continuous-time-Markov ML with a real rate matrix — so they wouldn't
|
||||
change how branch length is attributed per edge either.
|
||||
|
||||
**Export format note (not a bug in the generator).** Neither script's
|
||||
`.tre` output opens in PearTree (FigTree's successor) via File > Open —
|
||||
association/Launch-Services quirks were ruled out (the file was opened
|
||||
directly through the app, not by double-click). Likely cause, not yet
|
||||
confirmed: TNT's export is a minimal NEXUS `begin trees;` block with no
|
||||
preceding `Taxa` block and bare numeric (untranslated) leaf labels; PhyG's
|
||||
is multiple raw Newick trees concatenated with no NEXUS wrapper at all
|
||||
plus a trailing `[cost]` bracket tag after the root label. Both differ
|
||||
from a "normal" single-tree, fully-declared NEXUS file; this is PhyG/TNT's
|
||||
own export format, not something `write_sankoff_tnt`/`write_sankoff_phyg`
|
||||
could fix without post-processing the *other* program's output after the
|
||||
fact.
|
||||
|
||||
### Next direction: genuine ML branch lengths, not parsimony (open, 2026-08-12)
|
||||
|
||||
Decided: parsimony (the whole `--sankoff`/`--tnt`/`--phyg` pipeline above)
|
||||
is a stopgap, not the destination. The goal is maximum likelihood with
|
||||
real, calibrated branch lengths (expected substitutions/site), which
|
||||
parsimony step-counts were never going to give directly (see the open
|
||||
"branch lengths are `0.0`" problem above — even if fixed, TNT/PhyG-style
|
||||
parsimony branch length is a step count, not a continuous ML estimate).
|
||||
|
||||
**Model choices, settled:**
|
||||
- **The exchangeability `R` is symmetric; the rate matrix `Q` is not.**
|
||||
(Superseded an earlier, wrong framing here that treated "symmetric
|
||||
model" as one thing — see the resolution below on `R` vs `Q` vs `π` for
|
||||
the full reasoning.) `R(a,b) = R(b,a)` because `BasePairTally` never
|
||||
captured direction — a fact about the data, not a modelling choice.
|
||||
`Q(a,b) = R(a,b)·π_b` is asymmetric whenever the real state frequencies
|
||||
`π` are (which they are, empirically) — biology drives this via `π`,
|
||||
not via `R`.
|
||||
- **`∅` stays an ordinary 16th state**, as already established for
|
||||
TNT/PhyG — same reasoning applies to any ML tool: encode as a real
|
||||
alphabet symbol, never as `-`/gap, or the RAxML-era failure (empty set
|
||||
silently treated as missing data) repeats.
|
||||
|
||||
**Stationary frequencies for the 16 states — resolved (2026-08-12).**
|
||||
A CTMC needs a rate matrix `Q`, generally asymmetric. `Q(i,j) = R(i,j) ·
|
||||
π_j`, where `R` (exchangeability) is symmetric and `π` (stationary
|
||||
frequencies) need not be — this factoring is what makes `Q` reversible
|
||||
(satisfies detailed balance, `π_i·Q(i,j) = π_j·Q(j,i)`) for *any* `π`, not
|
||||
just uniform, as long as `R` is symmetric. Two separate, both-easy
|
||||
quantities, not one hard inverse problem:
|
||||
- **`R` is already calibrated**: `sub_cost` (`-ln(observed rate)` per pair,
|
||||
from `BasePairTally`) *is* `R` up to a log transform — recover it as
|
||||
`R(a,b) = exp(-sub_cost(a,b))`. Symmetric by construction, because the
|
||||
tally itself never distinguished direction (unordered-pair counts only)
|
||||
— not a modelling choice, a fact about what the data can say.
|
||||
- **`π` is a direct count**: empirical marginal frequency of each of the
|
||||
16 states across the whole alignment (same kind of scan already used to
|
||||
confirm `N` occurs 1383 times in the real 20-genome benchmark). With
|
||||
~908k sites × 20 genomes, the counts are large enough that this is
|
||||
precise on its own — no need to spend ML degrees of freedom
|
||||
re-estimating it via IQ-TREE's `+FO`.
|
||||
|
||||
Checked and ruled out along the way: IQ-TREE's `+F` (empirical, "compute
|
||||
from the alignment") does **not** work as a shortcut for this — for a
|
||||
custom-file morphology model, `readParameters` always requires the file's
|
||||
own frequency line unconditionally; omitting it and passing `+F` instead
|
||||
just fails (`ERROR: State frequencies could not be read`). `π` has to be
|
||||
computed by `obikmer` and written into the file, not left to IQ-TREE.
|
||||
|
||||
Net effect: no free-rate ML estimation needed for this piece at all (the
|
||||
mistaken assumption that motivated most of this discussion — that
|
||||
building an asymmetric `Q` a priori would require solving a linear system
|
||||
from `Q` itself — doesn't apply, because `R`, the only piece that's
|
||||
genuinely hard to get directionally, is symmetric and already in hand).
|
||||
|
||||
**Candidate tool: IQ-TREE**, because it supports user-defined multistate
|
||||
models (unlike RAxML's `MULTI` data type, which is limited to the
|
||||
equal-rate Mk model and can't take a custom rate matrix at all — a genuine
|
||||
tool limitation, not a gap-symbol encoding problem this time). IQ-TREE 3
|
||||
(3.0.1) is now installed locally (Homebrew, `iqtree3`).
|
||||
|
||||
**IQ-TREE custom-model format — verified empirically against the local
|
||||
binary (2026-08-12).** The web docs' `-mdef` NEXUS `begin models;
|
||||
frequency NAME = ...; model NAME = ...; end;` mechanism (initially assumed
|
||||
to apply directly, see history below) turned out to be for **named
|
||||
components used inside `MIX{...}`/`FMIX{...}` mixture models only** — it
|
||||
does **not** apply to a single, non-mixture custom morphology matrix, and
|
||||
using it that way fails (`ERROR: File not found <name>` — traced in
|
||||
IQ-TREE 3's own source, `model/modelmorphology.cpp`: any `-m` string that
|
||||
isn't `MK`/`ORDERED`/`GTR`/`GTRX` is passed straight to
|
||||
`ModelMarkov::readParameters()`, which opens it **as a literal file path**,
|
||||
never consulting the `-mdef` models block at all for this data type).
|
||||
|
||||
**The confirmed working recipe** (built a tiny 5-taxon/3-state toy dataset
|
||||
and rate file, ran it end to end with `iqtree3`, got a real ML tree with
|
||||
non-zero branch lengths and an optimized log-likelihood — ground truth,
|
||||
not documentation):
|
||||
- No `-mdef` needed. Write one plain file (any name) containing, as
|
||||
whitespace/newline-separated numbers, in order: the **lower-triangular
|
||||
rate matrix** (`N(N-1)/2` values, PAML row-major order — for 16 states,
|
||||
120 values, the same count and layout already produced for TNT's
|
||||
`smatrix`), immediately followed by the **N state frequencies** on the
|
||||
same stream (no header, no blank line required — confirmed by reading
|
||||
`ModelMorphology::readRates`/`ModelMarkov::readStateFreq` directly, which
|
||||
just pull tokens off the stream in sequence).
|
||||
- Invoke with `-m <path-to-that-file>+ASC` (`+ASC` for the no-constant-site
|
||||
correction, as before). An explicit `+F{f1,...,fN}` on the command line
|
||||
overrides the file's own frequency line if given (confirmed in
|
||||
`ModelMorphology::init`) — useful once real calibrated stationary
|
||||
frequencies exist, a placeholder equal-frequency line works meanwhile
|
||||
(the still-open gap noted above).
|
||||
- `--seqtype MORPH` (alphabet `0`-`9`,`A`-`V`, ≤31 states) — reuse the same
|
||||
`0-9A-F` recoding already built for TNT (`TNT_STATE_SYMBOL`).
|
||||
|
||||
**Risk, confirmed, precisely characterised, and resolved by design
|
||||
(2026-08-12).**
|
||||
`--seqtype MORPH{16}` **does not force the state count** for real ML
|
||||
analysis — tested directly (`--seqtype MORPH{4}` on the 3-symbol toy
|
||||
alignment gave the byte-for-byte identical 3-state result as no `{4}` at
|
||||
all) and confirmed in source: the value it sets
|
||||
(`params.alisim_num_states_morph`, `utils/tools.cpp`) is consumed only by
|
||||
the `--alisim` simulator; the main analysis path always calls
|
||||
`getDataBlockMorphStates`/an equivalent scan (`alignment/alignment.cpp`),
|
||||
for both FASTA/PHYLIP and NEXUS input (a NEXUS `symbols=` declaration
|
||||
doesn't change this either — checked, same code path). No CLI flag or
|
||||
NEXUS declaration overrides it.
|
||||
|
||||
The precise rule (from `getDataBlockMorphStates`, `alignment.cpp:1058`):
|
||||
`N` = **one plus the highest state ordinal actually observed anywhere in
|
||||
the alignment**, ordinal being the symbol's position in IQ-TREE's own
|
||||
fixed table `"0123456789ABCDEFGHIJKLMNOPQRSTUV"` — not a count of distinct
|
||||
symbols seen. So the risk is narrower than "any missing symbol breaks it":
|
||||
concretely, it's whether the symbol mapped to state index 15 (`F` in the
|
||||
`0-9A-F` recoding already used for TNT, i.e. our `N` = "all four bases
|
||||
ambiguous") occurs **at least once anywhere** in the real alignment — if
|
||||
it does, `N` correctly comes out to 16 regardless of which lower-index
|
||||
symbols (including `0`/`∅`) are rare or absent; if it doesn't, `N` silently
|
||||
undercounts and misaligns every value in a 16-entry rate/frequency file,
|
||||
with no error to catch it. **Checked against the real biological alignment**
|
||||
(`/tmp/msg_test/eub_sankoff.fasta`, the 20-genome benchmark index): the
|
||||
symbol `N` (IUPAC "all four bases ambiguous," state index 15) occurs 1383
|
||||
times across 13 of the 20 sequences — present, so this specific real
|
||||
dataset is not at risk. Still worth a real, general presence check inside
|
||||
`obikmer` before this is wired in, rather than assuming every future
|
||||
dataset will have `N` too (nothing in IQ-TREE would catch it if not).
|
||||
|
||||
**Resolution: subset + compact-renumber, not rely on all 16 appearing.**
|
||||
Since IQ-TREE always infers `N` from the alignment's own content and
|
||||
nothing overrides that, the fix is to make the file `obikmer` writes match
|
||||
that inference *by construction*, for every run, rather than hope the
|
||||
16th (or any particular) state happens to occur:
|
||||
1. Scan the real alignment for which of the 16 canonical states actually
|
||||
occur anywhere (not per-column — anywhere in the whole alignment).
|
||||
2. Renumber the occurring states to a **compact, consecutive** `0..k-1`
|
||||
range, preserving their relative order (the original bitmask/`STATE_
|
||||
SYMBOL` ordering) — not just filtering, since a *gap* in the ordinal
|
||||
sequence (e.g. keeping states `{0,1,2,4}` numbered as-is instead of
|
||||
`{0,1,2,3}`) reproduces the exact same "highest observed ordinal"
|
||||
miscount this was meant to fix.
|
||||
3. Recode the alignment itself with this new compact `k`-symbol alphabet
|
||||
(same recoding mechanism already used for TNT's `0-9A-F`, just over a
|
||||
possibly-smaller symbol set).
|
||||
4. Extract the matching `k×k` submatrix (rows/columns for the kept
|
||||
states only) from the full calibrated 16×16 cost matrix, in the same
|
||||
lower-triangular order the rate-matrix file needs — and, later, the
|
||||
matching `k`-length subset of stationary frequencies once those are
|
||||
calibrated (still the open gap noted earlier in this section).
|
||||
|
||||
Consequence, and why nothing is lost: a state that never occurs in a given
|
||||
alignment can, by definition, never contribute a transition to score in
|
||||
that same alignment — dropping it from that run's matrix costs nothing.
|
||||
The subset (and therefore `k`) can differ from one dataset/run to the
|
||||
next; this has to be done freshly per alignment, not computed once and
|
||||
reused.
|
||||
|
||||
*(Superseded reasoning, kept for the record: the `-mdef` NEXUS route
|
||||
below was the original plan, based on IQ-TREE's own web documentation for
|
||||
protein mixture models, before the local install allowed testing it —
|
||||
`GTRX` combined with a `-mdef`-referenced custom model, `+Fname` frequency
|
||||
reference. Both pieces exist and parse without error individually, but
|
||||
`GTRX`/`GTR` are IQ-TREE's own fixed built-in equal-structure multistate
|
||||
model, not a hook for an arbitrary custom matrix; a custom matrix is a
|
||||
file path in `-m` directly, no `-mdef` or `GTRX` involved.)*
|
||||
|
||||
Source for the empirical findings above: `model/modelmorphology.cpp` and
|
||||
`model/modelmarkov.cpp` in the local `iqtree/iqtree3` source (cloned to
|
||||
inspect the exact parsing logic after documentation didn't resolve the
|
||||
`+Fname` reference error) — more reliable here than the PDF/web manual,
|
||||
which (like TNT/PhyG) doesn't always match this specific binary. Original
|
||||
(partially superseded) sources: [Substitution
|
||||
Models](https://iqtree.github.io/doc/Substitution-Models), [Complex
|
||||
Models](https://iqtree.github.io/doc/Complex-Models).
|
||||
|
||||
**Relation to the existing calibration.** `sub_cost[a][b] = -ln(rate)`
|
||||
(see "A concrete Sankoff cost matrix" above) is already a log-rate — a
|
||||
genuine CTMC rate matrix `Q` could plausibly be recovered as
|
||||
`rate(a,b) = exp(-cost(a,b))`, renormalised so each row sums to zero, once
|
||||
the stationary-frequency gap above is closed. Not yet attempted.
|
||||
|
||||
## Heterozygosity, ploidy, and consensus-assembly inputs
|
||||
|
||||
A within-genome multiplicity signal (more than one of the 4 central forms
|
||||
|
||||
@@ -1,3 +1,7 @@
|
||||
mod phyg;
|
||||
mod sankoff;
|
||||
mod tnt;
|
||||
|
||||
use std::io::{self, BufWriter, Write};
|
||||
use std::path::PathBuf;
|
||||
|
||||
@@ -5,13 +9,17 @@ use clap::Args;
|
||||
use kodama::{Method, linkage};
|
||||
use obifastwrite::{JsonVal, write_record};
|
||||
use obikindex::{
|
||||
BasePairTally, DistanceMetric, KmerIndex, PHatEstimate, RawSnpDistanceOutput, SankoffWeights,
|
||||
DistanceMetric, KmerIndex, RawSnpDistanceOutput, SankoffWeights,
|
||||
SiblingAnnexStats, SnpAlignment, build_cost_matrix, c_ctx_from_p_hat, calibrate_p_hat,
|
||||
mean_substitution_cost, substitution_costs_from_tally,
|
||||
};
|
||||
use speedytree::{DistanceMatrix, Hybrid, NeighborJoiningSolver, to_newick};
|
||||
use tracing::info;
|
||||
|
||||
use phyg::write_sankoff_phyg;
|
||||
use sankoff::{write_sankoff_alignment_fasta, write_sankoff_matrix_csv, write_sankoff_params};
|
||||
use tnt::write_sankoff_tnt;
|
||||
|
||||
#[derive(clap::ValueEnum, Clone, Copy, Debug)]
|
||||
pub enum MetricArg {
|
||||
Jaccard,
|
||||
@@ -468,387 +476,6 @@ fn write_snp_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option<
|
||||
if n_sites == 1 { "" } else { "s" });
|
||||
}
|
||||
|
||||
// ── Sankoff pseudo-alignment → FASTA ────────────────────────────────────────
|
||||
//
|
||||
// Same data as `--snp`'s pseudo-alignment (`SnpAlignment`/
|
||||
// `snp_pseudo_alignment`), re-coded so its symbols match the accompanying
|
||||
// `--sankoff-matrix` output exactly: `0` for the empty/absent state instead
|
||||
// of `-`, which TNT/PhyG would otherwise read as their own gap character
|
||||
// rather than our "family absent" state.
|
||||
|
||||
fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option<PathBuf>) {
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff.fasta", p.display()))
|
||||
.unwrap_or_else(|| "sankoff.fasta".into());
|
||||
let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0);
|
||||
for (label, seq) in labels.iter().zip(alignment.sequences.iter()) {
|
||||
let recoded: Vec<u8> = seq.iter().map(|&b| if b == b'-' { b'0' } else { b }).collect();
|
||||
write_record(&recoded, label, &[("n_sites", JsonVal::Num(n_sites as u64))], &mut f).unwrap_or_else(|e| {
|
||||
eprintln!("error writing {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
}
|
||||
info!("Sankoff pseudo-alignment → {path} ({n_sites} site{})",
|
||||
if n_sites == 1 { "" } else { "s" });
|
||||
}
|
||||
|
||||
// ── Sankoff cost matrix → CSV ────────────────────────────────────────────────
|
||||
//
|
||||
// 16 states indexed by bitmask (bit 0=A, 1=C, 2=G, 3=T; state 0 is `∅`),
|
||||
// matching the convention already used for `--snp`'s IUPAC-coded output and
|
||||
// for the external TNT/PhyG scripts this feeds. Calibration report (p_hat,
|
||||
// its variance, how many pairs/loci went into it, the resulting c_ctx) goes
|
||||
// to the log, not the CSV, since it's a run-level fact, not per-cell data.
|
||||
|
||||
// IUPAC ambiguity code per state (same mapping as `siblings::iupac_code`,
|
||||
// already used for `--snp`'s pseudo-alignment — a biologist reads "R" as
|
||||
// "A or G" without needing this file's convention explained), with `0`
|
||||
// standing in for the empty state (`-` would collide with TNT/PhyG's own
|
||||
// gap/range syntax). Bit order: 0=A, 1=C, 2=G, 3=T. This project's
|
||||
// canonical alphabet — `write_sankoff_tnt` below recodes it to TNT's own
|
||||
// default alphabet at the adapter boundary, rather than using it here.
|
||||
const STATE_SYMBOL: [char; 16] = [
|
||||
'0', 'A', 'C', 'M', 'G', 'R', 'S', 'V', 'T', 'W', 'Y', 'H', 'K', 'D', 'B', 'N',
|
||||
];
|
||||
|
||||
fn write_sankoff_matrix_csv(
|
||||
matrix: &[[f64; 16]; 16],
|
||||
estimate: &PHatEstimate,
|
||||
weights: &SankoffWeights,
|
||||
output: &Option<PathBuf>,
|
||||
) {
|
||||
info!(
|
||||
p_hat = format_args!("{:.6}", estimate.p_hat),
|
||||
variance = format_args!("{:.3e}", estimate.variance),
|
||||
n_pairs_included = estimate.n_pairs_included,
|
||||
n_loci_total = estimate.n_loci_total,
|
||||
c_ctx = format_args!("{:.4}", weights.c_ctx),
|
||||
"Sankoff matrix calibration"
|
||||
);
|
||||
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff_matrix.csv", p.display()))
|
||||
.unwrap_or_else(|| "sankoff_matrix.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);
|
||||
}));
|
||||
write!(f, "state").unwrap();
|
||||
for sym in STATE_SYMBOL { write!(f, ",{sym}").unwrap(); }
|
||||
writeln!(f).unwrap();
|
||||
for (s, row) in matrix.iter().enumerate() {
|
||||
write!(f, "{}", STATE_SYMBOL[s]).unwrap();
|
||||
for cost in row { write!(f, ",{cost:.4}").unwrap(); }
|
||||
writeln!(f).unwrap();
|
||||
}
|
||||
info!("Sankoff cost matrix → {path}");
|
||||
}
|
||||
|
||||
// ── Sankoff calibration parameters → YAML report ────────────────────────────
|
||||
//
|
||||
// Everything `--sankoff` estimates from real data, in one durable,
|
||||
// machine-readable file: `p_hat` and its variance (with how many pairs/loci
|
||||
// went into it), the derived `c_ctx`, and the base-pair substitution tally
|
||||
// (raw counts, not just the derived costs) — kept for the same reason raw
|
||||
// counts are kept anywhere else in this project: costs are a modelling
|
||||
// choice built *from* the counts, and reproducing/re-deriving them later
|
||||
// needs the counts, not just their current derived value. Structured (YAML,
|
||||
// not an ad hoc key=value text file) so R/Python/etc. can load it directly
|
||||
// rather than re-parsing free text.
|
||||
|
||||
#[derive(serde::Serialize)]
|
||||
struct SankoffSubstitution {
|
||||
pair: String,
|
||||
count: u64,
|
||||
cost: f64,
|
||||
}
|
||||
|
||||
#[derive(serde::Serialize)]
|
||||
struct SankoffParamsReport {
|
||||
ratio_ceiling: f64,
|
||||
flank_length_m: usize,
|
||||
p_hat: f64,
|
||||
p_hat_variance: f64,
|
||||
n_pairs_included: usize,
|
||||
n_loci_total: u64,
|
||||
/// Weighted-average substitution cost — see `c_ctx_from_p_hat`'s docs:
|
||||
/// this is what turns the raw expected mutation *count* behind `c_ctx`
|
||||
/// into an actual cost (not every mutation is worth a flat `1`).
|
||||
mean_sub_cost: f64,
|
||||
c_ctx: f64,
|
||||
substitutions: Vec<SankoffSubstitution>,
|
||||
}
|
||||
|
||||
fn write_sankoff_params(
|
||||
estimate: &PHatEstimate,
|
||||
tally: &BasePairTally,
|
||||
weights: &SankoffWeights,
|
||||
ratio_ceiling: f64,
|
||||
m: usize,
|
||||
mean_sub_cost: f64,
|
||||
output: &Option<PathBuf>,
|
||||
) {
|
||||
const BASE_LETTER: [char; 4] = ['A', 'C', 'G', 'T'];
|
||||
|
||||
let mut substitutions = Vec::with_capacity(6);
|
||||
for a in 0..4 {
|
||||
for b in (a + 1)..4 {
|
||||
substitutions.push(SankoffSubstitution {
|
||||
pair: format!("{}/{}", BASE_LETTER[a], BASE_LETTER[b]),
|
||||
count: tally.counts[a][b],
|
||||
cost: weights.sub_cost[a][b],
|
||||
});
|
||||
}
|
||||
}
|
||||
let report = SankoffParamsReport {
|
||||
ratio_ceiling,
|
||||
flank_length_m: m,
|
||||
p_hat: estimate.p_hat,
|
||||
p_hat_variance: estimate.variance,
|
||||
n_pairs_included: estimate.n_pairs_included,
|
||||
n_loci_total: estimate.n_loci_total,
|
||||
mean_sub_cost,
|
||||
c_ctx: weights.c_ctx,
|
||||
substitutions,
|
||||
};
|
||||
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff_params.yaml", p.display()))
|
||||
.unwrap_or_else(|| "sankoff_params.yaml".into());
|
||||
let f = std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
serde_yaml::to_writer(f, &report).unwrap_or_else(|e| {
|
||||
eprintln!("error writing {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
info!("Sankoff calibration parameters → {path}");
|
||||
}
|
||||
|
||||
// ── Sankoff cost matrix + alignment → ready-to-run TNT script ──────────────
|
||||
//
|
||||
// TNT's *default* xread reader only accepts its own 0-9A-F alphabet (see
|
||||
// its manual: "up to 16 states are allowed by xread, using symbols 0-9 ...
|
||||
// and A-F") — the wider IUPAC set `STATE_SYMBOL` uses is rejected as an
|
||||
// "alien symbol" unless `nstates dna` is set, which imposes TNT's own fixed
|
||||
// DNA encoding instead, incompatible with a custom smatrix. And TNT's
|
||||
// `smatrix`/`cost` commands reject decimal costs ("found symbol . when
|
||||
// reading transformation costs"). So: recode to TNT's alphabet and
|
||||
// integer-scale the costs here, at this adapter's boundary, rather than
|
||||
// degrading the project's own canonical (IUPAC, real-valued) output.
|
||||
|
||||
const TNT_STATE_SYMBOL: [char; 16] = [
|
||||
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F',
|
||||
];
|
||||
|
||||
fn write_sankoff_tnt(
|
||||
matrix: &[[f64; 16]; 16],
|
||||
alignment: &SnpAlignment,
|
||||
labels: &[String],
|
||||
output: &Option<PathBuf>,
|
||||
cost_scale: f64,
|
||||
) {
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff.tnt", p.display()))
|
||||
.unwrap_or_else(|| "sankoff.tnt".into());
|
||||
let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
|
||||
// IUPAC-ish symbol -> bitmask, to translate the alignment (which uses
|
||||
// `STATE_SYMBOL`, `-` already normalised to `0` by `snp_pseudo_alignment`
|
||||
// callers) into TNT's alphabet without re-deriving state indices.
|
||||
let mut iupac_to_state = [0u8; 128];
|
||||
for (state, &sym) in STATE_SYMBOL.iter().enumerate() {
|
||||
iupac_to_state[sym as usize] = state as u8;
|
||||
}
|
||||
|
||||
let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0);
|
||||
writeln!(f, "xread").unwrap();
|
||||
writeln!(f, "'obikmer central-position SNP families, calibrated Sankoff 16-state encoding'").unwrap();
|
||||
writeln!(f, "{n_sites} {}", labels.len()).unwrap();
|
||||
for (label, seq) in labels.iter().zip(alignment.sequences.iter()) {
|
||||
write!(f, "{label} ").unwrap();
|
||||
for &b in seq {
|
||||
let b = if b == b'-' { b'0' } else { b };
|
||||
let state = iupac_to_state[b as usize];
|
||||
write!(f, "{}", TNT_STATE_SYMBOL[state as usize]).unwrap();
|
||||
}
|
||||
writeln!(f).unwrap();
|
||||
}
|
||||
writeln!(f, ";\n").unwrap();
|
||||
|
||||
let scaled_matrix = scaled_metric_matrix(matrix, cost_scale);
|
||||
|
||||
writeln!(f, "smatrix =0 (family16)").unwrap();
|
||||
for i in 0..16 {
|
||||
for j in (i + 1)..16 {
|
||||
writeln!(f, "{}/{} {}", TNT_STATE_SYMBOL[i], TNT_STATE_SYMBOL[j], scaled_matrix[i][j]).unwrap();
|
||||
}
|
||||
}
|
||||
writeln!(f, ";\n").unwrap();
|
||||
|
||||
writeln!(f, "ccode ( 0.{} ;", n_sites - 1).unwrap();
|
||||
writeln!(f, "smatrix +0 0.{} ;", n_sites - 1).unwrap();
|
||||
writeln!(f).unwrap();
|
||||
|
||||
// Basename only (not the full `path`/`output` prefix): TNT's natural
|
||||
// workflow is to `cd` into the output directory before `proc`-ing the
|
||||
// script, and an absolute path here would break if that directory is
|
||||
// later moved or copied elsewhere.
|
||||
let tre_name = output.as_ref()
|
||||
.and_then(|p| p.file_name())
|
||||
.map(|n| format!("{}_sankoff.tre", n.to_string_lossy()))
|
||||
.unwrap_or_else(|| "sankoff.tre".into());
|
||||
|
||||
// TNT's plain command parser has no comment syntax of its own — `/* */`
|
||||
// and `[ ]` are only recognised inside the (separately-enabled) macro
|
||||
// scripting language, and fail with "No command!" here otherwise
|
||||
// (verified against this file with the local TNT binary). `quote` is
|
||||
// the closest working equivalent: it prints free text and does not
|
||||
// otherwise affect parsing, so it doubles as an explanation of the
|
||||
// defaults below when the script is run. `;` ends a `quote` block like
|
||||
// any other TNT command, so the text itself must avoid semicolons.
|
||||
writeln!(f, "quote").unwrap();
|
||||
writeln!(f, "Default search below (edit or delete this block to run your own strategy):").unwrap();
|
||||
writeln!(f, " hold N : size of TNT's tree buffer (how many equally-parsimonious").unwrap();
|
||||
writeln!(f, " trees it keeps in memory at once), 20 is a small, fast").unwrap();
|
||||
writeln!(f, " default, raise it if mult reports it had to drop trees.").unwrap();
|
||||
writeln!(f, " mult : traditional search (random addition sequences followed by").unwrap();
|
||||
writeln!(f, " TBR branch-swapping, TNT's own default replication count),").unwrap();
|
||||
writeln!(f, " a reasonable first-pass strategy on this data's memory").unwrap();
|
||||
writeln!(f, " footprint, xmult's ratchet/drift/tree-fusion buffers ran").unwrap();
|
||||
writeln!(f, " this out of RAM at TNT's default mxram on this dataset.").unwrap();
|
||||
writeln!(f, " export - F : write the trees held in the buffer to file F, in").unwrap();
|
||||
writeln!(f, " TNT/Hennig86 format ('-' means trees, as opposed to data).").unwrap();
|
||||
writeln!(f, ";").unwrap();
|
||||
writeln!(f, "hold 20;").unwrap();
|
||||
writeln!(f, "mult;").unwrap();
|
||||
writeln!(f, "export - {tre_name};").unwrap();
|
||||
|
||||
info!(
|
||||
"TNT script → {path} (costs scaled x{cost_scale:.0}, runs a default `hold 20; mult;` \
|
||||
search and writes trees to {tre_name} in TNT's working directory — edit the trailing \
|
||||
comment block in the script to change this)\n\
|
||||
Run it with:\n \
|
||||
printf 'proc {path};\\nquit;\\n' | tnt\n\
|
||||
(or start `tnt` interactively and type `proc {path};`)"
|
||||
);
|
||||
}
|
||||
|
||||
/// Scale `matrix` by `cost_scale` and round to integers (TNT's smatrix/cost
|
||||
/// commands reject decimals), then take the *metric closure* of the result
|
||||
/// (Floyd-Warshall over the 16 states again, on the now-integer values).
|
||||
///
|
||||
/// `matrix` is already a metric in its real-valued form (it's a
|
||||
/// shortest-path closure itself — see `build_cost_matrix`), but rounding
|
||||
/// each cell independently can still break the triangle inequality: e.g.
|
||||
/// two real costs of `1.734` each round to `173`, summing to `346`, while
|
||||
/// their own real sum `3.468` rounds to `347` — TNT then reports "triangle
|
||||
/// inequality violated ... Fixed" and silently substitutes its own
|
||||
/// corrected value. Re-closing after rounding makes that correction
|
||||
/// explicit and reproducible here instead, rather than left implicit and
|
||||
/// TNT-version-dependent.
|
||||
fn scaled_metric_matrix(matrix: &[[f64; 16]; 16], cost_scale: f64) -> [[i64; 16]; 16] {
|
||||
let mut m = [[0i64; 16]; 16];
|
||||
for i in 0..16 {
|
||||
for j in 0..16 {
|
||||
m[i][j] = (matrix[i][j] * cost_scale).round() as i64;
|
||||
}
|
||||
}
|
||||
for k in 0..16 {
|
||||
for i in 0..16 {
|
||||
for j in 0..16 {
|
||||
let via = m[i][k] + m[k][j];
|
||||
if via < m[i][j] {
|
||||
m[i][j] = via;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
// ── Sankoff cost matrix → PhyG custom-alphabet TCM + ready-to-run script ────
|
||||
//
|
||||
// PhyG's `tcm:STRING` format needs no alphabet recoding, unlike `--tnt`:
|
||||
// its parser reads the alphabet straight from the tcm file's own first
|
||||
// line, so `--sankoff`'s own `_sankoff.fasta` (already IUPAC+`0`) is reused
|
||||
// as-is via `prefasta:`. PhyG auto-adds its own indel/gap state as an
|
||||
// (n+1)-th row/column of the tcm — inert here since the alignment already
|
||||
// encodes absence as an ordinary state (`0`), never as `-` (see
|
||||
// `write_sankoff_alignment_fasta`'s own comment on why, and the RAxML-era
|
||||
// bug that motivated it). The gap row/column below reuses `matrix[i][0]`/
|
||||
// `matrix[0][j]` (cost to/from `∅`) as the closest principled value for a
|
||||
// state that, in practice, is never actually triggered.
|
||||
|
||||
fn write_sankoff_phyg(matrix: &[[f64; 16]; 16], output: &Option<PathBuf>, cost_scale: f64) {
|
||||
let scaled_matrix = scaled_metric_matrix(matrix, cost_scale);
|
||||
|
||||
let basename = |suffix: &str| -> String {
|
||||
output.as_ref()
|
||||
.and_then(|p| p.file_name())
|
||||
.map(|n| format!("{}{suffix}", n.to_string_lossy()))
|
||||
.unwrap_or_else(|| format!("sankoff{suffix}"))
|
||||
};
|
||||
let full_path = |suffix: &str| -> String {
|
||||
output.as_ref()
|
||||
.map(|p| format!("{}{suffix}", p.display()))
|
||||
.unwrap_or_else(|| format!("sankoff{suffix}"))
|
||||
};
|
||||
|
||||
let tcm_path = full_path("_sankoff.tcm");
|
||||
let mut f = BufWriter::new(std::fs::File::create(&tcm_path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {tcm_path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let alphabet_line = STATE_SYMBOL.iter().map(|c| c.to_string()).collect::<Vec<_>>().join(" ");
|
||||
writeln!(f, "{alphabet_line}").unwrap();
|
||||
for i in 0..16 {
|
||||
let mut row: Vec<i64> = (0..16).map(|j| scaled_matrix[i][j]).collect();
|
||||
row.push(scaled_matrix[i][0]); // gap column: same cost as to/from ∅
|
||||
writeln!(f, "{}", row.iter().map(|v| v.to_string()).collect::<Vec<_>>().join(" ")).unwrap();
|
||||
}
|
||||
let mut gap_row: Vec<i64> = (0..16).map(|j| scaled_matrix[0][j]).collect();
|
||||
gap_row.push(0);
|
||||
writeln!(f, "{}", gap_row.iter().map(|v| v.to_string()).collect::<Vec<_>>().join(" ")).unwrap();
|
||||
info!("PhyG TCM → {tcm_path}");
|
||||
|
||||
let pg_path = full_path("_sankoff.pg");
|
||||
let mut f = BufWriter::new(std::fs::File::create(&pg_path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {pg_path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let fasta_name = basename("_sankoff.fasta");
|
||||
let tcm_name = basename("_sankoff.tcm");
|
||||
let tre_name = basename("_sankoff.tre");
|
||||
writeln!(f, "read(prefasta:\"{fasta_name}\", tcm:\"{tcm_name}\")").unwrap();
|
||||
writeln!(f, "search(seconds:300, instances:4)").unwrap();
|
||||
writeln!(f, "report(\"{tre_name}\", graphs, newick, overwrite)").unwrap();
|
||||
|
||||
let pg_dir = std::path::Path::new(&pg_path).parent()
|
||||
.filter(|d| !d.as_os_str().is_empty())
|
||||
.map(|d| d.display().to_string())
|
||||
.unwrap_or_else(|| ".".into());
|
||||
let pg_name = std::path::Path::new(&pg_path).file_name()
|
||||
.map(|n| n.to_string_lossy().into_owned())
|
||||
.unwrap_or_else(|| pg_path.clone());
|
||||
info!(
|
||||
"PhyG script → {pg_path} (costs scaled x{cost_scale:.0}, runs a default 300s/4-instance \
|
||||
search and writes trees to {tre_name})\n\
|
||||
Run it with:\n \
|
||||
cd {pg_dir} && phyg {pg_name}\n\
|
||||
(`phyg` must run from that directory — `read()`/`report()` in the script use relative \
|
||||
file names)"
|
||||
);
|
||||
}
|
||||
|
||||
// ── UPGMA Newick from kodama dendrogram ───────────────────────────────────────
|
||||
|
||||
fn upgma_to_newick(dendro: &kodama::Dendrogram<f64>, names: &[String]) -> String {
|
||||
@@ -0,0 +1,81 @@
|
||||
use std::io::{BufWriter, Write};
|
||||
use std::path::PathBuf;
|
||||
|
||||
use tracing::info;
|
||||
|
||||
use super::sankoff::{STATE_SYMBOL, scaled_metric_matrix};
|
||||
|
||||
// ── Sankoff cost matrix → PhyG custom-alphabet TCM + ready-to-run script ────
|
||||
//
|
||||
// PhyG's `tcm:STRING` format needs no alphabet recoding, unlike `--tnt`:
|
||||
// its parser reads the alphabet straight from the tcm file's own first
|
||||
// line, so `--sankoff`'s own `_sankoff.fasta` (already IUPAC+`0`) is reused
|
||||
// as-is via `prefasta:`. PhyG auto-adds its own indel/gap state as an
|
||||
// (n+1)-th row/column of the tcm — inert here since the alignment already
|
||||
// encodes absence as an ordinary state (`0`), never as `-` (see
|
||||
// `sankoff::write_sankoff_alignment_fasta`'s own comment on why, and the
|
||||
// RAxML-era bug that motivated it). The gap row/column below reuses
|
||||
// `matrix[i][0]`/`matrix[0][j]` (cost to/from `∅`) as the closest
|
||||
// principled value for a state that, in practice, is never actually
|
||||
// triggered.
|
||||
|
||||
pub(super) fn write_sankoff_phyg(matrix: &[[f64; 16]; 16], output: &Option<PathBuf>, cost_scale: f64) {
|
||||
let scaled_matrix = scaled_metric_matrix(matrix, cost_scale);
|
||||
|
||||
let basename = |suffix: &str| -> String {
|
||||
output.as_ref()
|
||||
.and_then(|p| p.file_name())
|
||||
.map(|n| format!("{}{suffix}", n.to_string_lossy()))
|
||||
.unwrap_or_else(|| format!("sankoff{suffix}"))
|
||||
};
|
||||
let full_path = |suffix: &str| -> String {
|
||||
output.as_ref()
|
||||
.map(|p| format!("{}{suffix}", p.display()))
|
||||
.unwrap_or_else(|| format!("sankoff{suffix}"))
|
||||
};
|
||||
|
||||
let tcm_path = full_path("_sankoff.tcm");
|
||||
let mut f = BufWriter::new(std::fs::File::create(&tcm_path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {tcm_path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let alphabet_line = STATE_SYMBOL.iter().map(|c| c.to_string()).collect::<Vec<_>>().join(" ");
|
||||
writeln!(f, "{alphabet_line}").unwrap();
|
||||
for i in 0..16 {
|
||||
let mut row: Vec<i64> = (0..16).map(|j| scaled_matrix[i][j]).collect();
|
||||
row.push(scaled_matrix[i][0]); // gap column: same cost as to/from ∅
|
||||
writeln!(f, "{}", row.iter().map(|v| v.to_string()).collect::<Vec<_>>().join(" ")).unwrap();
|
||||
}
|
||||
let mut gap_row: Vec<i64> = (0..16).map(|j| scaled_matrix[0][j]).collect();
|
||||
gap_row.push(0);
|
||||
writeln!(f, "{}", gap_row.iter().map(|v| v.to_string()).collect::<Vec<_>>().join(" ")).unwrap();
|
||||
info!("PhyG TCM → {tcm_path}");
|
||||
|
||||
let pg_path = full_path("_sankoff.pg");
|
||||
let mut f = BufWriter::new(std::fs::File::create(&pg_path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {pg_path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let fasta_name = basename("_sankoff.fasta");
|
||||
let tcm_name = basename("_sankoff.tcm");
|
||||
let tre_name = basename("_sankoff.tre");
|
||||
writeln!(f, "read(prefasta:\"{fasta_name}\", tcm:\"{tcm_name}\")").unwrap();
|
||||
writeln!(f, "search(seconds:300, instances:4)").unwrap();
|
||||
writeln!(f, "report(\"{tre_name}\", graphs, newick, overwrite)").unwrap();
|
||||
|
||||
let pg_dir = std::path::Path::new(&pg_path).parent()
|
||||
.filter(|d| !d.as_os_str().is_empty())
|
||||
.map(|d| d.display().to_string())
|
||||
.unwrap_or_else(|| ".".into());
|
||||
let pg_name = std::path::Path::new(&pg_path).file_name()
|
||||
.map(|n| n.to_string_lossy().into_owned())
|
||||
.unwrap_or_else(|| pg_path.clone());
|
||||
info!(
|
||||
"PhyG script → {pg_path} (costs scaled x{cost_scale:.0}, runs a default 300s/4-instance \
|
||||
search and writes trees to {tre_name})\n\
|
||||
Run it with:\n \
|
||||
cd {pg_dir} && phyg {pg_name}\n\
|
||||
(`phyg` must run from that directory — `read()`/`report()` in the script use relative \
|
||||
file names)"
|
||||
);
|
||||
}
|
||||
@@ -0,0 +1,202 @@
|
||||
use std::io::{BufWriter, Write};
|
||||
use std::path::PathBuf;
|
||||
|
||||
use obifastwrite::{JsonVal, write_record};
|
||||
use obikindex::{BasePairTally, PHatEstimate, SankoffWeights, SnpAlignment};
|
||||
use tracing::info;
|
||||
|
||||
// ── Sankoff pseudo-alignment → FASTA ────────────────────────────────────────
|
||||
//
|
||||
// Same data as `--snp`'s pseudo-alignment (`SnpAlignment`/
|
||||
// `snp_pseudo_alignment`), re-coded so its symbols match the accompanying
|
||||
// `--sankoff-matrix` output exactly: `0` for the empty/absent state instead
|
||||
// of `-`, which TNT/PhyG would otherwise read as their own gap character
|
||||
// rather than our "family absent" state.
|
||||
|
||||
pub(super) fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option<PathBuf>) {
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff.fasta", p.display()))
|
||||
.unwrap_or_else(|| "sankoff.fasta".into());
|
||||
let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0);
|
||||
for (label, seq) in labels.iter().zip(alignment.sequences.iter()) {
|
||||
let recoded: Vec<u8> = seq.iter().map(|&b| if b == b'-' { b'0' } else { b }).collect();
|
||||
write_record(&recoded, label, &[("n_sites", JsonVal::Num(n_sites as u64))], &mut f).unwrap_or_else(|e| {
|
||||
eprintln!("error writing {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
}
|
||||
info!("Sankoff pseudo-alignment → {path} ({n_sites} site{})",
|
||||
if n_sites == 1 { "" } else { "s" });
|
||||
}
|
||||
|
||||
// ── Sankoff cost matrix → CSV ────────────────────────────────────────────────
|
||||
//
|
||||
// 16 states indexed by bitmask (bit 0=A, 1=C, 2=G, 3=T; state 0 is `∅`),
|
||||
// matching the convention already used for `--snp`'s IUPAC-coded output and
|
||||
// for the external TNT/PhyG scripts this feeds. Calibration report (p_hat,
|
||||
// its variance, how many pairs/loci went into it, the resulting c_ctx) goes
|
||||
// to the log, not the CSV, since it's a run-level fact, not per-cell data.
|
||||
|
||||
// IUPAC ambiguity code per state (same mapping as `siblings::iupac_code`,
|
||||
// already used for `--snp`'s pseudo-alignment — a biologist reads "R" as
|
||||
// "A or G" without needing this file's convention explained), with `0`
|
||||
// standing in for the empty state (`-` would collide with TNT/PhyG's own
|
||||
// gap/range syntax). Bit order: 0=A, 1=C, 2=G, 3=T. This project's
|
||||
// canonical alphabet — `tnt::write_sankoff_tnt` recodes it to TNT's own
|
||||
// default alphabet at the adapter boundary, rather than using it here.
|
||||
pub(super) const STATE_SYMBOL: [char; 16] = [
|
||||
'0', 'A', 'C', 'M', 'G', 'R', 'S', 'V', 'T', 'W', 'Y', 'H', 'K', 'D', 'B', 'N',
|
||||
];
|
||||
|
||||
pub(super) fn write_sankoff_matrix_csv(
|
||||
matrix: &[[f64; 16]; 16],
|
||||
estimate: &PHatEstimate,
|
||||
weights: &SankoffWeights,
|
||||
output: &Option<PathBuf>,
|
||||
) {
|
||||
info!(
|
||||
p_hat = format_args!("{:.6}", estimate.p_hat),
|
||||
variance = format_args!("{:.3e}", estimate.variance),
|
||||
n_pairs_included = estimate.n_pairs_included,
|
||||
n_loci_total = estimate.n_loci_total,
|
||||
c_ctx = format_args!("{:.4}", weights.c_ctx),
|
||||
"Sankoff matrix calibration"
|
||||
);
|
||||
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff_matrix.csv", p.display()))
|
||||
.unwrap_or_else(|| "sankoff_matrix.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);
|
||||
}));
|
||||
write!(f, "state").unwrap();
|
||||
for sym in STATE_SYMBOL { write!(f, ",{sym}").unwrap(); }
|
||||
writeln!(f).unwrap();
|
||||
for (s, row) in matrix.iter().enumerate() {
|
||||
write!(f, "{}", STATE_SYMBOL[s]).unwrap();
|
||||
for cost in row { write!(f, ",{cost:.4}").unwrap(); }
|
||||
writeln!(f).unwrap();
|
||||
}
|
||||
info!("Sankoff cost matrix → {path}");
|
||||
}
|
||||
|
||||
// ── Sankoff calibration parameters → YAML report ────────────────────────────
|
||||
//
|
||||
// Everything `--sankoff` estimates from real data, in one durable,
|
||||
// machine-readable file: `p_hat` and its variance (with how many pairs/loci
|
||||
// went into it), the derived `c_ctx`, and the base-pair substitution tally
|
||||
// (raw counts, not just the derived costs) — kept for the same reason raw
|
||||
// counts are kept anywhere else in this project: costs are a modelling
|
||||
// choice built *from* the counts, and reproducing/re-deriving them later
|
||||
// needs the counts, not just their current derived value. Structured (YAML,
|
||||
// not an ad hoc key=value text file) so R/Python/etc. can load it directly
|
||||
// rather than re-parsing free text.
|
||||
|
||||
#[derive(serde::Serialize)]
|
||||
struct SankoffSubstitution {
|
||||
pair: String,
|
||||
count: u64,
|
||||
cost: f64,
|
||||
}
|
||||
|
||||
#[derive(serde::Serialize)]
|
||||
struct SankoffParamsReport {
|
||||
ratio_ceiling: f64,
|
||||
flank_length_m: usize,
|
||||
p_hat: f64,
|
||||
p_hat_variance: f64,
|
||||
n_pairs_included: usize,
|
||||
n_loci_total: u64,
|
||||
/// Weighted-average substitution cost — see `c_ctx_from_p_hat`'s docs:
|
||||
/// this is what turns the raw expected mutation *count* behind `c_ctx`
|
||||
/// into an actual cost (not every mutation is worth a flat `1`).
|
||||
mean_sub_cost: f64,
|
||||
c_ctx: f64,
|
||||
substitutions: Vec<SankoffSubstitution>,
|
||||
}
|
||||
|
||||
pub(super) fn write_sankoff_params(
|
||||
estimate: &PHatEstimate,
|
||||
tally: &BasePairTally,
|
||||
weights: &SankoffWeights,
|
||||
ratio_ceiling: f64,
|
||||
m: usize,
|
||||
mean_sub_cost: f64,
|
||||
output: &Option<PathBuf>,
|
||||
) {
|
||||
const BASE_LETTER: [char; 4] = ['A', 'C', 'G', 'T'];
|
||||
|
||||
let mut substitutions = Vec::with_capacity(6);
|
||||
for a in 0..4 {
|
||||
for b in (a + 1)..4 {
|
||||
substitutions.push(SankoffSubstitution {
|
||||
pair: format!("{}/{}", BASE_LETTER[a], BASE_LETTER[b]),
|
||||
count: tally.counts[a][b],
|
||||
cost: weights.sub_cost[a][b],
|
||||
});
|
||||
}
|
||||
}
|
||||
let report = SankoffParamsReport {
|
||||
ratio_ceiling,
|
||||
flank_length_m: m,
|
||||
p_hat: estimate.p_hat,
|
||||
p_hat_variance: estimate.variance,
|
||||
n_pairs_included: estimate.n_pairs_included,
|
||||
n_loci_total: estimate.n_loci_total,
|
||||
mean_sub_cost,
|
||||
c_ctx: weights.c_ctx,
|
||||
substitutions,
|
||||
};
|
||||
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff_params.yaml", p.display()))
|
||||
.unwrap_or_else(|| "sankoff_params.yaml".into());
|
||||
let f = std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
serde_yaml::to_writer(f, &report).unwrap_or_else(|e| {
|
||||
eprintln!("error writing {path}: {e}");
|
||||
std::process::exit(1);
|
||||
});
|
||||
info!("Sankoff calibration parameters → {path}");
|
||||
}
|
||||
|
||||
/// Scale `matrix` by `cost_scale` and round to integers (TNT's smatrix/cost
|
||||
/// and PhyG's `tcm:` commands both reject decimals), then take the *metric
|
||||
/// closure* of the result (Floyd-Warshall over the 16 states again, on the
|
||||
/// now-integer values).
|
||||
///
|
||||
/// `matrix` is already a metric in its real-valued form (it's a
|
||||
/// shortest-path closure itself — see `obikindex::build_cost_matrix`), but
|
||||
/// rounding each cell independently can still break the triangle
|
||||
/// inequality: e.g. two real costs of `1.734` each round to `173`, summing
|
||||
/// to `346`, while their own real sum `3.468` rounds to `347` — TNT then
|
||||
/// reports "triangle inequality violated ... Fixed" and silently
|
||||
/// substitutes its own corrected value. Re-closing after rounding makes
|
||||
/// that correction explicit and reproducible here instead, rather than
|
||||
/// left implicit and tool-version-dependent.
|
||||
pub(super) fn scaled_metric_matrix(matrix: &[[f64; 16]; 16], cost_scale: f64) -> [[i64; 16]; 16] {
|
||||
let mut m = [[0i64; 16]; 16];
|
||||
for i in 0..16 {
|
||||
for j in 0..16 {
|
||||
m[i][j] = (matrix[i][j] * cost_scale).round() as i64;
|
||||
}
|
||||
}
|
||||
for k in 0..16 {
|
||||
for i in 0..16 {
|
||||
for j in 0..16 {
|
||||
let via = m[i][k] + m[k][j];
|
||||
if via < m[i][j] {
|
||||
m[i][j] = via;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
@@ -0,0 +1,119 @@
|
||||
use std::io::{BufWriter, Write};
|
||||
use std::path::PathBuf;
|
||||
|
||||
use obikindex::SnpAlignment;
|
||||
use tracing::info;
|
||||
|
||||
use super::sankoff::{STATE_SYMBOL, scaled_metric_matrix};
|
||||
|
||||
// ── Sankoff cost matrix + alignment → ready-to-run TNT script ──────────────
|
||||
//
|
||||
// TNT's *default* xread reader only accepts its own 0-9A-F alphabet (see
|
||||
// its manual: "up to 16 states are allowed by xread, using symbols 0-9 ...
|
||||
// and A-F") — the wider IUPAC set `STATE_SYMBOL` uses is rejected as an
|
||||
// "alien symbol" unless `nstates dna` is set, which imposes TNT's own fixed
|
||||
// DNA encoding instead, incompatible with a custom smatrix. And TNT's
|
||||
// `smatrix`/`cost` commands reject decimal costs ("found symbol . when
|
||||
// reading transformation costs"). So: recode to TNT's alphabet and
|
||||
// integer-scale the costs here, at this adapter's boundary, rather than
|
||||
// degrading the project's own canonical (IUPAC, real-valued) output.
|
||||
|
||||
const TNT_STATE_SYMBOL: [char; 16] = [
|
||||
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F',
|
||||
];
|
||||
|
||||
pub(super) fn write_sankoff_tnt(
|
||||
matrix: &[[f64; 16]; 16],
|
||||
alignment: &SnpAlignment,
|
||||
labels: &[String],
|
||||
output: &Option<PathBuf>,
|
||||
cost_scale: f64,
|
||||
) {
|
||||
let path = output.as_ref()
|
||||
.map(|p| format!("{}_sankoff.tnt", p.display()))
|
||||
.unwrap_or_else(|| "sankoff.tnt".into());
|
||||
let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| {
|
||||
eprintln!("error creating {path}: {e}");
|
||||
std::process::exit(1);
|
||||
}));
|
||||
|
||||
// IUPAC-ish symbol -> bitmask, to translate the alignment (which uses
|
||||
// `STATE_SYMBOL`, `-` already normalised to `0` by `snp_pseudo_alignment`
|
||||
// callers) into TNT's alphabet without re-deriving state indices.
|
||||
let mut iupac_to_state = [0u8; 128];
|
||||
for (state, &sym) in STATE_SYMBOL.iter().enumerate() {
|
||||
iupac_to_state[sym as usize] = state as u8;
|
||||
}
|
||||
|
||||
let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0);
|
||||
writeln!(f, "xread").unwrap();
|
||||
writeln!(f, "'obikmer central-position SNP families, calibrated Sankoff 16-state encoding'").unwrap();
|
||||
writeln!(f, "{n_sites} {}", labels.len()).unwrap();
|
||||
for (label, seq) in labels.iter().zip(alignment.sequences.iter()) {
|
||||
write!(f, "{label} ").unwrap();
|
||||
for &b in seq {
|
||||
let b = if b == b'-' { b'0' } else { b };
|
||||
let state = iupac_to_state[b as usize];
|
||||
write!(f, "{}", TNT_STATE_SYMBOL[state as usize]).unwrap();
|
||||
}
|
||||
writeln!(f).unwrap();
|
||||
}
|
||||
writeln!(f, ";\n").unwrap();
|
||||
|
||||
let scaled_matrix = scaled_metric_matrix(matrix, cost_scale);
|
||||
|
||||
writeln!(f, "smatrix =0 (family16)").unwrap();
|
||||
for i in 0..16 {
|
||||
for j in (i + 1)..16 {
|
||||
writeln!(f, "{}/{} {}", TNT_STATE_SYMBOL[i], TNT_STATE_SYMBOL[j], scaled_matrix[i][j]).unwrap();
|
||||
}
|
||||
}
|
||||
writeln!(f, ";\n").unwrap();
|
||||
|
||||
writeln!(f, "ccode ( 0.{} ;", n_sites - 1).unwrap();
|
||||
writeln!(f, "smatrix +0 0.{} ;", n_sites - 1).unwrap();
|
||||
writeln!(f).unwrap();
|
||||
|
||||
// Basename only (not the full `path`/`output` prefix): TNT's natural
|
||||
// workflow is to `cd` into the output directory before `proc`-ing the
|
||||
// script, and an absolute path here would break if that directory is
|
||||
// later moved or copied elsewhere.
|
||||
let tre_name = output.as_ref()
|
||||
.and_then(|p| p.file_name())
|
||||
.map(|n| format!("{}_sankoff.tre", n.to_string_lossy()))
|
||||
.unwrap_or_else(|| "sankoff.tre".into());
|
||||
|
||||
// TNT's plain command parser has no comment syntax of its own — `/* */`
|
||||
// and `[ ]` are only recognised inside the (separately-enabled) macro
|
||||
// scripting language, and fail with "No command!" here otherwise
|
||||
// (verified against this file with the local TNT binary). `quote` is
|
||||
// the closest working equivalent: it prints free text and does not
|
||||
// otherwise affect parsing, so it doubles as an explanation of the
|
||||
// defaults below when the script is run. `;` ends a `quote` block like
|
||||
// any other TNT command, so the text itself must avoid semicolons.
|
||||
writeln!(f, "quote").unwrap();
|
||||
writeln!(f, "Default search below (edit or delete this block to run your own strategy):").unwrap();
|
||||
writeln!(f, " hold N : size of TNT's tree buffer (how many equally-parsimonious").unwrap();
|
||||
writeln!(f, " trees it keeps in memory at once), 20 is a small, fast").unwrap();
|
||||
writeln!(f, " default, raise it if mult reports it had to drop trees.").unwrap();
|
||||
writeln!(f, " mult : traditional search (random addition sequences followed by").unwrap();
|
||||
writeln!(f, " TBR branch-swapping, TNT's own default replication count),").unwrap();
|
||||
writeln!(f, " a reasonable first-pass strategy on this data's memory").unwrap();
|
||||
writeln!(f, " footprint, xmult's ratchet/drift/tree-fusion buffers ran").unwrap();
|
||||
writeln!(f, " this out of RAM at TNT's default mxram on this dataset.").unwrap();
|
||||
writeln!(f, " export - F : write the trees held in the buffer to file F, in").unwrap();
|
||||
writeln!(f, " TNT/Hennig86 format ('-' means trees, as opposed to data).").unwrap();
|
||||
writeln!(f, ";").unwrap();
|
||||
writeln!(f, "hold 20;").unwrap();
|
||||
writeln!(f, "mult;").unwrap();
|
||||
writeln!(f, "export - {tre_name};").unwrap();
|
||||
|
||||
info!(
|
||||
"TNT script → {path} (costs scaled x{cost_scale:.0}, runs a default `hold 20; mult;` \
|
||||
search and writes trees to {tre_name} in TNT's working directory — edit the trailing \
|
||||
comment block in the script to change this)\n\
|
||||
Run it with:\n \
|
||||
printf 'proc {path};\\nquit;\\n' | tnt\n\
|
||||
(or start `tnt` interactively and type `proc {path};`)"
|
||||
);
|
||||
}
|
||||
Reference in New Issue
Block a user