diff --git a/docmd/theory/evolutionary_distances.md b/docmd/theory/evolutionary_distances.md index b1fd239f..6fd21d8b 100644 --- a/docmd/theory/evolutionary_distances.md +++ b/docmd/theory/evolutionary_distances.md @@ -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 ` — 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 +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 diff --git a/src/obikmer/src/cmd/distance.rs b/src/obikmer/src/cmd/distance/mod.rs similarity index 55% rename from src/obikmer/src/cmd/distance.rs rename to src/obikmer/src/cmd/distance/mod.rs index bd7a4e37..bd698caa 100644 --- a/src/obikmer/src/cmd/distance.rs +++ b/src/obikmer/src/cmd/distance/mod.rs @@ -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) { - 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 = 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, -) { - 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, -} - -fn write_sankoff_params( - estimate: &PHatEstimate, - tally: &BasePairTally, - weights: &SankoffWeights, - ratio_ceiling: f64, - m: usize, - mean_sub_cost: f64, - output: &Option, -) { - 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, - 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, 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::>().join(" "); - writeln!(f, "{alphabet_line}").unwrap(); - for i in 0..16 { - let mut row: Vec = (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::>().join(" ")).unwrap(); - } - let mut gap_row: Vec = (0..16).map(|j| scaled_matrix[0][j]).collect(); - gap_row.push(0); - writeln!(f, "{}", gap_row.iter().map(|v| v.to_string()).collect::>().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, names: &[String]) -> String { diff --git a/src/obikmer/src/cmd/distance/phyg.rs b/src/obikmer/src/cmd/distance/phyg.rs new file mode 100644 index 00000000..f0d0167d --- /dev/null +++ b/src/obikmer/src/cmd/distance/phyg.rs @@ -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, 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::>().join(" "); + writeln!(f, "{alphabet_line}").unwrap(); + for i in 0..16 { + let mut row: Vec = (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::>().join(" ")).unwrap(); + } + let mut gap_row: Vec = (0..16).map(|j| scaled_matrix[0][j]).collect(); + gap_row.push(0); + writeln!(f, "{}", gap_row.iter().map(|v| v.to_string()).collect::>().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)" + ); +} diff --git a/src/obikmer/src/cmd/distance/sankoff.rs b/src/obikmer/src/cmd/distance/sankoff.rs new file mode 100644 index 00000000..76676237 --- /dev/null +++ b/src/obikmer/src/cmd/distance/sankoff.rs @@ -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) { + 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 = 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, +) { + 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, +} + +pub(super) fn write_sankoff_params( + estimate: &PHatEstimate, + tally: &BasePairTally, + weights: &SankoffWeights, + ratio_ceiling: f64, + m: usize, + mean_sub_cost: f64, + output: &Option, +) { + 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 +} diff --git a/src/obikmer/src/cmd/distance/tnt.rs b/src/obikmer/src/cmd/distance/tnt.rs new file mode 100644 index 00000000..977994c4 --- /dev/null +++ b/src/obikmer/src/cmd/distance/tnt.rs @@ -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, + 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};`)" + ); +}