From 8615da59a84a035846321e529b02e2163d96a38a Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Thu, 13 Aug 2026 21:14:35 +0200 Subject: [PATCH] Add phylogenetic CLI options for family overlap and missing data Introduces CLI flags for computing pairwise family overlap matrices and filtering genomes below a shared family threshold. Adds a free-loss mode that recodes locus non-detection states to missing data symbols in Sankoff-calibrated alignments, resolving ascertainment bias handling for IQ-TREE. Updates empirical transition parameters, removes the legacy model asset, and extends output writers for CSV diagnostics, FASTA pseudo-alignments, and Newick trees. --- .gitignore | 9 + UserDocMD/usage/phylo.md | 13 +- docmd/theory/evolutionary_distances.md | 120 ++++++ iqtree.model | 16 - sankoff_params.yaml | 148 ++++---- src/obikmer/src/cmd/phylo/args.rs | 243 ++++++++++++ src/obikmer/src/cmd/phylo/family_overlap.rs | 97 +++++ src/obikmer/src/cmd/phylo/iqtree.rs | 73 +++- src/obikmer/src/cmd/phylo/mod.rs | 390 ++------------------ src/obikmer/src/cmd/phylo/outputs.rs | 168 +++++++++ src/obikmer/src/cmd/phylo/sankoff.rs | 13 +- src/obikmer/src/cmd/phylo/tnt.rs | 8 + 12 files changed, 833 insertions(+), 465 deletions(-) delete mode 100644 iqtree.model create mode 100644 src/obikmer/src/cmd/phylo/args.rs create mode 100644 src/obikmer/src/cmd/phylo/family_overlap.rs create mode 100644 src/obikmer/src/cmd/phylo/outputs.rs diff --git a/.gitignore b/.gitignore index 32189eb2..addd10f0 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,12 @@ benchmark/specific_index_presence TNT phyg *.tnt +*.tre +*.phy +*.treefile +*.bionj +*.iqtree +*.mldist +*.parstree +*.ckp.gz +*.model diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md index 4a71bd0c..265c5767 100644 --- a/UserDocMD/usage/phylo.md +++ b/UserDocMD/usage/phylo.md @@ -47,10 +47,12 @@ A **family** is the set of up to 4 kmers that share identical flanking sequence |---|---|---| | `--sibling-annex` | off | Build the sibling presence-mask annex (prerequisite for every option below) | | `--exclude-genome LABEL` | none | Exclude a genome (repeatable) from every SNP/Sankoff/export computation below | +| `--min-shared-family N` | none | Auto-exclude any genome whose mean shared-family count against every other genome (see `--family-overlap`) falls below `N` — same exclusion as `--exclude-genome`, applied on top of it | | `--sibling-stats` | off | Write the family-size (sibling count) distribution, per genome and globally | | `--raw-snp-distance` | off | Write the single-copy central-SNP p-distance matrix | | `--raw-snp-counts` | off | Write per-pair diagnostic counts (n_snp, n_shared, n_eligible) instead of a matrix | | `--snp` | off | Write a SNP-only pseudo-alignment in FASTA, IUPAC-coded | +| `--family-overlap` | off | Write an NxN matrix of, for each genome pair, how many variable families both genomes actually carry a call for; the diagonal holds each genome's own total family count | ### Locus eligibility @@ -62,12 +64,19 @@ A family is eligible for a genome pair $(i, j)$ only if genome $i$ carries exact `--exclude-genome` removes a genome from these computations, re-checking column variability among the remaining genomes so that a column made monomorphic by the exclusion is dropped rather than kept artificially. It does not affect the `--metric` distance-matrix path. +### Family overlap and low-coverage genomes + +`--family-overlap` writes, for every genome pair, how many variable families both genomes actually carry a call for (neither is absent) — a direct measure of how much informative content two genomes actually share. On genome-skim or otherwise incomplete-coverage collections, a genome with very little overlap with everything else has almost nothing left to constrain its position in a tree, and tends to end up placed unstably (near-zero branch length, grafted inside an unrelated clade) by `--tnt`/`--iqtree`. + +`--min-shared-family N` automates the fix: it excludes, before any computation, every genome whose mean shared-family count against all other genomes (the same statistic, averaged per row of the `--family-overlap` matrix) falls below `N`. There is no universal value for `N` — it depends on how divergent and how completely covered the genome collection is; inspect `--family-overlap`'s own output to find where the real gap sits before choosing a threshold. + ## Sankoff calibration and phylogenetic exports | Option | Default | Description | |---|---|---| | `--sankoff` | off | Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment | | `--sankoff-ratio-ceiling` | `0.5` | Exclude genome pairs whose raw SNP ratio exceeds this value from the calibration | +| `--free-loss` | off | Recode a family's non-detection as the `?` missing-data symbol instead of an ordinary, costed state, in `--sankoff`'s pseudo-alignment and every export built from it | | `--tnt` | off | Also write a TNT script (implies `--sankoff`) | | `--phyg` | off | Also write PhyG input files (implies `--sankoff`) | | `--iqtree` | off | Also write an IQ-TREE custom model and alignment (implies `--sankoff`) | @@ -84,6 +93,8 @@ which are combined into a row-normalized $16 \times 16$ transition probability m `--sankoff` alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol `0` (never a gap character, to avoid ambiguity with external tools' own gap semantics). It does not run any external tool. +With `--free-loss`, the empty state is recoded to `?` instead — TNT/PhyG/IQ-TREE's own missing-data symbol — rather than an ordinary, costed 16th state. This matters for genome-skim or otherwise incomplete-coverage collections, where non-detection of a family is dominated by sampling failure rather than true evolutionary loss: scoring it as a real state risks grouping genomes by shared undersampling instead of shared ancestry. `?` rather than `-` because `-` still carries gap/indel semantics in these tools, and a non-detected family is not an observed deletion. Combine with `--min-shared-family`/`--family-overlap` above: `--free-loss` removes the false signal from non-detection, but a genome left with too little real overlap with everything else will still be placed unstably — excluding it is the other half of the fix. + ### Exports All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment, recoded for the target tool: @@ -94,4 +105,4 @@ All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment, ## Output files -With `-o/--output PREFIX`, the relevant subset of the following files is written: `_dist.csv`, `_shared.csv`, `_nj.nwk`, `_upgma.nwk`, `_siblings.csv`, `_rawsnp.csv`, `_rawsnp_counts.csv`, `_snp.fasta`, `_sankoff_matrix.csv`, `_sankoff_params.yaml`, `_sankoff.fasta`, `_sankoff.tnt`, `_sankoff.tcm`, `_sankoff.pg`, `_iqtree.model`, `_iqtree.fasta`. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout. +With `-o/--output PREFIX`, the relevant subset of the following files is written: `_dist.csv`, `_shared.csv`, `_nj.nwk`, `_upgma.nwk`, `_siblings.csv`, `_rawsnp.csv`, `_rawsnp_counts.csv`, `_snp.fasta`, `_family_overlap.csv`, `_sankoff_matrix.csv`, `_sankoff_params.yaml`, `_sankoff.fasta`, `_sankoff.tnt`, `_sankoff.tcm`, `_sankoff.pg`, `_iqtree.model`, `_iqtree.fasta`. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout. diff --git a/docmd/theory/evolutionary_distances.md b/docmd/theory/evolutionary_distances.md index bf5a3be6..1bb5d1a9 100644 --- a/docmd/theory/evolutionary_distances.md +++ b/docmd/theory/evolutionary_distances.md @@ -1836,6 +1836,126 @@ lookups proportionally; `p_hat` is unbiased. Off by default (exact). matrix). 5. Subsampling (Step 5). +## Locus dropout under incomplete coverage: absence as missing data, not a penalized state (2026-08-13, discussion) + +Distinct from the SNP/Sankoff calibration above (which conditions on +*conserved-flank* comparability), this addresses the coarser presence/absence +signal of a family itself, for genome-skim/reduced-representation inputs +where per-genome coverage is often `<1x`, sometimes `<0.5x`. At that depth, +non-detection of a locus is dominated by Poisson sampling failure, not by +true loss — even for a single-copy locus, `P(non-detection | present)` is +non-negligible. Treating "absent" as a real 0-state therefore risks +clustering genomes by shared undersampling rather than shared ancestry. + +**Parsimony.** Loss and gain are the two directions of the same transition, +so they cannot be freed asymmetrically without a biological justification for +trusting one direction and not the other — none exists here. Freeing both +directions to zero cost makes the character contribute nothing to any tree's +score regardless of topology, which is operationally equivalent to just +coding the state as ambiguous (`?`) rather than 0/1. `?` is the more honest +encoding when downstream ancestral-state reconstruction is wanted: cost-zero +Fitch/Sankoff still assigns internal states arbitrarily among equal-cost +options, whereas `?` makes the "we don't know" explicit. + +**ML.** The correct analogue is not a tuned asymmetric loss/gain rate but +literal missing-data coding at the tip: an ambiguous state (all-states tip +likelihood vector, e.g. `(1,1)` for a 2-state character) so Felsenstein +pruning marginalises over it and the character contributes likelihood 1 +(neutral) for that taxon — the same mechanism already used for the `N` +IUPAC-ambiguous case in the 16-state Sankoff encoding (see "IQ-TREE +custom-model format" above). No rate parameter to justify or estimate. + +**Open item, not yet decided.** No per-locus/per-genome confidence +distinction is implemented: the current plan is to code *every* non-detected +locus as missing/ambiguous, uniformly, without a depth-based threshold for +"confidently absent." Deferred because at typical skim coverage the +confidently-absent case is rare enough that a uniform treatment is +defensible as a first pass. Known consequence to watch for: a genome with +very low overall coverage ends up with a mostly-`?` row and constrains the +tree weakly — will likely need a minimum-detected-loci-per-genome filter +upstream (exclude or flag under-covered genomes) before this becomes a +correctness problem rather than a power problem. Not designed yet. + +### `--free-loss` implemented (2026-08-13) + +`obikmer phylo --sankoff --free-loss` (with `--tnt`/`--phyg`/`--iqtree`) +recodes `∅` (a family's non-detection in a genome) to `?` instead of the +default ordinary, costed 16th state — the resolution of the discussion +above. `?`, not `-`: `-` still carries gap/indel semantics in TNT/PhyG/ +IQ-TREE (see "`∅` stays an ordinary 16th state" above, and the reasoning +that motivated it in the first place — the RAxML-era failure where `-` +was silently swallowed as missing data was undesired back when `∅` was +meant to be fully informative; `--free-loss` is the first case where that +same tool behaviour is exactly what's wanted). `?` is the standard +missing/ambiguous symbol these tools already recognise without needing to +be declared in the custom alphabet, so parsimony treats it as +"compatible with any state" (Sankoff/Fitch's own missing-data extension — +no per-transition cost edit needed) and Felsenstein pruning marginalises +it for IQ-TREE's ML. + +Implementation is a recoding at each export's existing `-`-normalisation +boundary, not a change to matrix calibration: `write_sankoff_alignment_fasta` +(`sankoff.rs`, feeds `--phyg` too via its `prefasta:`), `write_sankoff_tnt` +(`tnt.rs`, `?` written directly, bypassing `TNT_STATE_SYMBOL` — there is no +state index for it), and `compact_alphabet`/`write_iqtree_alignment` +(`iqtree.rs`, `∅` excluded from the compact alphabet entirely rather than +assigned a symbol). The calibrated cost matrix itself is unchanged — its +`∅`-row/column simply becomes unused once no sequence contains that state. + +**Unverified, flagged rather than assumed** (this project's own standard for +tool-specific claims — see the TNT/PhyG/IQ-TREE syntax facts above, all +checked against the real binaries, not the manuals alone): +- PhyG auto-adds its own `(n+1)`-th gap/indel row to the `tcm:` matrix + (`phyg.rs`); with `--free-loss` the reused `_sankoff.fasta` now contains + real `?` characters PhyG did not have to interpret before. Not confirmed + against the local PhyG binary that `?` is read as "any state" rather than + routed through that gap row, or otherwise mishandled. +- IQ-TREE's exact `?`/missing-data semantics for `--seqtype MORPH` with a + custom-file model plus `+ASC` — assumed to marginalise via Felsenstein + pruning by the general convention, not yet re-verified against the local + `iqtree3` build the way the rest of the custom-model mechanism was. + +Both should be checked against a real run (small index) before `--free-loss` +results are trusted for TNT/PhyG/IQ-TREE, the same way every other +tool-specific behaviour in this file was established. + +### TNT confirmed (2026-08-13); IQ-TREE `+ASC` constant-column fix (2026-08-13) + +**TNT.** Real run, 20-genome benchmark index, `--free-loss --tnt`: `?` +present for all 20 taxa in both `sankoff.fasta` and `sankoff.tnt`, loaded +with no syntax error and no triangle-inequality warning, `mult` search +completed, `sankoff.tre` holds 20 equally-parsimonious trees. Topology is +exactly what the theory predicts: strain-level clades unaffected (the 4 *E. +coli*, the *Klebsiella* trio, the *Salmonella* pairs — driven by real +substitution signal, not gene-content) but resolution among the deeply +divergent lineages (*Candidozyma*, *Saccharolobus*, *Wolbachia*, +*Acidobacterium*, *Opitutus*, *Bacillus*/*Shouchella*) is weak and unstable +across the 20 tied trees — the expected trade-off of discarding the +gene-content/dropout signal that used to help place exactly those lineages. +TNT's own handling of `?` — read as "compatible with any state" — no longer +a flagged assumption. + +**IQ-TREE: `+ASC` broke on real data — non-informative columns survive +`--free-loss`, distinct bug from the `--exclude-genome` one above.** +`snp_pseudo_alignment`'s "variable family" test (`family_size() >= 2`) is a +whole-annex property, blind to any single column's actual calls. Once +non-detection is hidden as `?` (ignored for constancy) rather than counted +as its own real state, a family that is genuinely variable *somewhere in the +dataset* can still have only one distinct call among the genomes that +actually carry it in one particular column, with the rest `?` — constant +once missing data is excluded, which is exactly what `+ASC` forbids. Same +underlying failure mode as `--exclude-genome`'s (documented above, "Two +consistency bugs found and fixed post-implementation"), same remedy: rescan +columns and drop the ones that are constant once `-` is ignored — but a +different trigger (hiding cells vs. dropping whole rows), so a separate fix +was needed rather than reusing `drop_excluded` directly. Implemented as +`iqtree::drop_ascertainment_noninformative`, run only under `--free-loss`, +only on IQ-TREE's own copy of the alignment (`--tnt`/`--phyg` have no +no-invariant-site requirement and must not lose those columns). Not yet +re-run against the local `iqtree3` binary to confirm `+ASC` now accepts the +filtered alignment — next verification step, alongside the still-open PhyG +`?`-vs-gap-row question above. + ## References The Mash mutation-rate model this discussion contrasts with: diff --git a/iqtree.model b/iqtree.model deleted file mode 100644 index 09bd6ac5..00000000 --- a/iqtree.model +++ /dev/null @@ -1,16 +0,0 @@ -0.222641 -0.222347 0.005197 -0.206041 0.017442 0.019145 -0.222349 0.024185 0.003581 0.013346 -0.208691 0.017667 0.002921 0.000012 0.019168 -0.205784 0.013974 0.019121 0.000085 0.018901 0.000018 -0.027849 0.001121 0.001230 0.000109 0.001216 0.000109 0.000118 -0.222640 0.004745 0.023685 0.013831 0.005020 0.002986 0.013813 0.000889 -0.206337 0.017467 0.012967 0.000076 0.013365 0.000016 0.000059 0.000079 0.017590 -0.208662 0.003028 0.019388 0.000017 0.002764 0.000003 0.000018 0.000019 0.017788 0.000017 -0.027866 0.001122 0.001231 0.000109 0.000858 0.000016 0.000083 0.000009 0.001129 0.000100 0.000111 -0.206081 0.013994 0.012951 0.000059 0.018929 0.000017 0.000082 0.000085 0.017568 0.000078 0.000012 0.000076 -0.027872 0.001122 0.000833 0.000074 0.001217 0.000109 0.000080 0.000042 0.001130 0.000100 0.000017 0.000006 0.000108 -0.027846 0.000899 0.001230 0.000087 0.001216 0.000019 0.000118 0.000009 0.001129 0.000080 0.000111 0.000044 0.000108 0.000009 -0.009724 0.000172 0.000189 0.000044 0.000187 0.000044 0.000048 0.000047 0.000174 0.000041 0.000045 0.000043 0.000044 0.000043 0.000047 -0.904811 0.012806 0.013642 0.007640 0.013786 0.007933 0.004196 0.000140 0.012382 0.006685 0.007847 0.000195 0.007527 0.000202 0.000134 0.000076 diff --git a/sankoff_params.yaml b/sankoff_params.yaml index e4a0075c..489cf29c 100644 --- a/sankoff_params.yaml +++ b/sankoff_params.yaml @@ -2,166 +2,166 @@ ratio_ceiling: 0.5 cardinality_transitions: - from: 0 to: 0 - count: 117158166 - probability: 0.8249649657257814 + count: 37622567 + probability: 0.7765461217831444 - from: 0 to: 1 - count: 13707223 - probability: 0.09651891232567299 + count: 6977400 + probability: 0.14401656617768033 - from: 0 to: 2 - count: 10934100 - probability: 0.07699206755884405 + count: 3744558 + probability: 0.0772893033240408 - from: 0 to: 3 - count: 193315 - probability: 0.0013612205430842902 + count: 91083 + probability: 0.0018799926759482985 - from: 0 to: 4 - count: 23125 - probability: 0.00016283384661730445 + count: 12985 + probability: 0.0002680160391861122 - from: 1 to: 0 - count: 13707223 - probability: 0.8980473700324103 + count: 6977400 + probability: 0.820630384894398 - from: 1 to: 1 - count: 1043692 - probability: 0.06837890181868832 + count: 1018042 + probability: 0.1197346000370715 - from: 1 to: 2 - count: 507953 - probability: 0.033279232106318904 + count: 502624 + probability: 0.05911493200578466 - from: 1 to: 3 - count: 4270 - probability: 0.0002797548613631216 + count: 4197 + probability: 0.000493620220340211 - from: 1 to: 4 count: 225 - probability: 0.000014741181219368235 + probability: 0.000026462842405658202 - from: 2 to: 0 - count: 10934100 - probability: 0.9551033885158036 + count: 3744558 + probability: 0.8804244781723005 - from: 2 to: 1 - count: 507953 - probability: 0.04437014765794788 + count: 502624 + probability: 0.11817749195415704 - from: 2 to: 2 - count: 5241 - probability: 0.00045780602511512846 + count: 5162 + probability: 0.0012136949558066439 - from: 2 to: 3 - count: 690 - probability: 0.000060272115498843476 + count: 688 + probability: 0.0001617632951559417 - from: 2 to: 4 count: 96 - probability: 8.385685634621701e-6 + probability: 0.00002257162257989884 - from: 3 to: 0 - count: 193315 - probability: 0.9743748708410829 + count: 91083 + probability: 0.9478728718311618 - from: 3 to: 1 - count: 4270 - probability: 0.021522285898618442 + count: 4197 + probability: 0.04367689297756317 - from: 3 to: 2 - count: 690 - probability: 0.0034778401100812 + count: 688 + probability: 0.007159805186696083 - from: 3 to: 3 count: 98 - probability: 0.0004939541025912429 + probability: 0.0010198559713607792 - from: 3 to: 4 count: 26 - probability: 0.0001310490476262481 + probability: 0.00027057403321816594 - from: 4 to: 0 - count: 23125 - probability: 0.9846291407647109 + count: 12985 + probability: 0.972950696838004 - from: 4 to: 1 count: 225 - probability: 0.009580175423656646 + probability: 0.01685898396523303 - from: 4 to: 2 count: 96 - probability: 0.004087541514093502 + probability: 0.0071931664918327585 - from: 4 to: 3 count: 26 - probability: 0.0011070424934003236 + probability: 0.0019481492582047056 - from: 4 to: 4 count: 14 - probability: 0.0005960998041386358 + probability: 0.0010490034467256107 composition_transitions: - from: 'A' to: 'A' - count: 162260 - probability: 0.47062167539692207 + count: 157182 + probability: 0.46648741052031767 - from: 'A' to: 'C' - count: 27810 - probability: 0.08066059899413536 + count: 27305 + probability: 0.08103624298111281 - from: 'A' to: 'G' - count: 130153 - probability: 0.3774979842101294 + count: 128548 + probability: 0.38150693875612857 - from: 'A' to: 'T' - count: 24555 - probability: 0.07121974139881315 + count: 23913 + probability: 0.07096940774244097 - from: 'C' to: 'A' - count: 27810 - probability: 0.07790832534920075 + count: 27305 + probability: 0.07823849488677553 - from: 'C' to: 'C' - count: 184633 - probability: 0.5172401234879174 + count: 179118 + probability: 0.5132365034656459 - from: 'C' to: 'G' - count: 19637 - probability: 0.05501207424963161 + count: 19392 + probability: 0.055564947549692406 - from: 'C' to: 'T' - count: 124878 - probability: 0.3498394769132503 + count: 123182 + probability: 0.3529600540978862 - from: 'G' to: 'A' - count: 130153 - probability: 0.36057158212891627 + count: 128548 + probability: 0.3637269077214298 - from: 'G' to: 'C' - count: 19637 - probability: 0.054401697680925745 + count: 19392 + probability: 0.05486971555009776 - from: 'G' to: 'G' - count: 184557 - probability: 0.5112906308956874 + count: 179245 + probability: 0.5071742039901647 - from: 'G' to: 'T' - count: 26616 - probability: 0.07373608929447062 + count: 26234 + probability: 0.07422917273830779 - from: 'T' to: 'A' - count: 24555 - probability: 0.07337692220342934 + count: 23913 + probability: 0.07307212790143376 - from: 'T' to: 'C' - count: 124878 - probability: 0.3731689387464813 + count: 123182 + probability: 0.3764132839524281 - from: 'T' to: 'G' - count: 26616 - probability: 0.07953574267426085 + count: 26234 + probability: 0.08016452153080807 - from: 'T' to: 'T' - count: 158593 - probability: 0.4739183963758285 + count: 153923 + probability: 0.4703500666153301 diff --git a/src/obikmer/src/cmd/phylo/args.rs b/src/obikmer/src/cmd/phylo/args.rs new file mode 100644 index 00000000..deedb5bd --- /dev/null +++ b/src/obikmer/src/cmd/phylo/args.rs @@ -0,0 +1,243 @@ +use std::path::PathBuf; + +use clap::Args; +use obikindex::DistanceMetric; + +#[derive(clap::ValueEnum, Clone, Copy, Debug)] +pub enum MetricArg { + Jaccard, + Mash, + Hamming, + BrayCurtis, + #[value(name = "relfreq-bray-curtis")] + RelfreqBrayCurtis, + Euclidean, + #[value(name = "relfreq-euclidean")] + RelfreqEuclidean, + Hellinger, + #[value(name = "hellinger-euclidean")] + HellingerEuclidean, +} + +impl From for DistanceMetric { + fn from(m: MetricArg) -> Self { + match m { + MetricArg::Jaccard => DistanceMetric::Jaccard, + MetricArg::Mash => DistanceMetric::Mash, + MetricArg::Hamming => DistanceMetric::Hamming, + MetricArg::BrayCurtis => DistanceMetric::BrayCurtis, + MetricArg::RelfreqBrayCurtis => DistanceMetric::RelfreqBrayCurtis, + MetricArg::Euclidean => DistanceMetric::Euclidean, + MetricArg::RelfreqEuclidean => DistanceMetric::RelfreqEuclidean, + MetricArg::Hellinger => DistanceMetric::Hellinger, + MetricArg::HellingerEuclidean => DistanceMetric::HellingerEuclidean, + } + } +} + +#[derive(Args)] +pub struct PhyloArgs { + /// Index directory + pub index: PathBuf, + + /// Distance metric to compute + #[arg(long, value_enum, default_value = "jaccard")] + pub metric: MetricArg, + + /// Minimum count to consider a kmer present when computing Jaccard on count indexes + #[arg(long, default_value = "1")] + pub presence_threshold: u32, + + /// Also output the shared-kmer count matrix (CSV) + #[arg(long)] + pub shared_kmers: bool, + + /// Compute and write a Neighbor-Joining tree (Newick) + #[arg(long)] + pub nj: bool, + + /// Compute and write a UPGMA tree (Newick) + #[arg(long)] + pub upgma: bool, + + /// Build the sibling-count/minorant annex on this (multi-genome) index + /// — see `docmd/theory/evolutionary_distances.md`, Step 2b. Construction + /// only; does not by itself compute or write any statistics. + #[arg(long)] + pub sibling_annex: bool, + + /// Exclude a genome (by its exact label) from every computation below + /// that reads the sibling annex — `--raw-snp-distance`/`--raw-snp-counts`, + /// `--snp`, and `--sankoff` (and everything `--sankoff` implies: the + /// cardinality/composition transition models, the exported + /// matrix/alignment, `--tnt`/`--phyg`/`--iqtree`). Repeatable. Does + /// *not* affect the plain `--metric` distance matrix/NJ/UPGMA path (a + /// different, unrelated computation). Applied by zeroing the excluded + /// genome's row/column after `raw_snp_distance` runs (a pair with zero + /// counts is already skipped by `base_pair_tally`/`cardinality_tally`, + /// so this needs no change to the underlying traversal) and by + /// dropping its row from `snp_pseudo_alignment`'s output — the annex + /// is still built/scanned for the excluded genome too, just not used + /// afterward. For a genome with almost no informative sites shared + /// with anything else (see `docmd/theory/evolutionary_distances.md`, + /// the IQ-TREE/Mash rogue-taxon discussion), its presence can + /// otherwise silently bias the transition models. + #[arg(long = "exclude-genome", value_name = "LABEL")] + pub exclude_genome: Vec, + + /// Auto-exclude any genome whose mean shared-family count against every + /// other genome (same statistic as `--family-overlap`'s matrix, averaged + /// over each row excluding the diagonal) falls below this threshold — + /// same exclusion machinery as `--exclude-genome`, applied on top of it + /// rather than instead of it. Empirically, genomes below ~1000 shared + /// families on the 20-genome benchmark are exactly the ones that placed + /// themselves arbitrarily under `--tnt`/`--iqtree` (near-zero branch + /// lengths, grafted inside unrelated clades) — too little real + /// constraint on where they belong. See + /// `docmd/theory/evolutionary_distances.md`, "Locus dropout under + /// incomplete coverage". + #[arg(long, value_name = "N")] + pub min_shared_family: Option, + + /// Tally the sibling-count distribution (CSV) of an already-built annex + /// (run with `--sibling-annex` first, in this invocation or an earlier + /// one). A separate, occasional diagnostic pass — not run every time the + /// annex itself is (re)built. + #[arg(long)] + pub sibling_stats: bool, + + /// Compute the raw p-distance restricted to loci that are single-copy + /// in both genomes of each pair (an already-built sibling annex is + /// required — run with `--sibling-annex` first, in this invocation or + /// an earlier one). A quick way to test the central-position SNP + /// estimator against a real index; not the full `SnpTally` design. + #[arg(long)] + pub raw_snp_distance: bool, + + /// Write the raw per-pair counts (`n_snp`, `n_shared`, `n_eligible`) + /// behind `--raw-snp-distance`'s ratio, one row per genome pair — a + /// diagnostic table, not a matrix. The ratio alone can't distinguish + /// "identical at every eligible locus" from "almost no eligible loci + /// at all" (e.g. `0.0` from 0/2 looks the same as `0.0` from 0/2000), + /// and that distinction matters a lot for genome pairs near the edge + /// of what central-position families can resolve (see + /// `docmd/theory/evolutionary_distances.md`, "Run 3" and the + /// IQ-TREE/Mash comparison). Same annex requirement as + /// `--raw-snp-distance`. + #[arg(long)] + pub raw_snp_counts: bool, + + /// Write a SNP-only pseudo-alignment (FASTA, IUPAC-coded) from an + /// already-built sibling annex — one row per genome, one column per + /// variable family (monomorphic families skipped), no flanking + /// sequence. See `docmd/theory/evolutionary_distances.md`, + /// "Multi-genome framing: family as pseudo-alignment column". + #[arg(long)] + pub snp: bool, + + /// Write an NxN CSV (`_family_overlap.csv`) of, for each genome + /// pair, how many variable families (same set `--snp`'s pseudo-alignment + /// uses — `family_size() >= 2`) both genomes actually carry a call for + /// (neither is `∅`). A direct read of how much informative content two + /// genomes actually share at the family level — the diagnostic for why + /// a genome with little overlap with anything else (e.g. an + /// under-covered or very divergent one) ends up placed unstably by + /// `--tnt`/`--iqtree`: little-to-no shared, real data to constrain it. + /// Same annex requirement as `--snp`. + #[arg(long)] + pub family_overlap: bool, + + /// Calibrate a 16-state Sankoff cost matrix and its matching + /// pseudo-alignment from an already-built sibling annex (run with + /// `--sibling-annex` first, in this invocation or an earlier one), for + /// use with TNT/PhyG. See `docmd/theory/evolutionary_distances.md`, + /// "Sankoff parsimony as the resolution of the 16-state model problem". + #[arg(long)] + pub sankoff: bool, + + /// Recode a family's non-detection (`∅`, no member observed in a + /// genome) as TNT/PhyG/IQ-TREE's own missing-data symbol (`?`) in + /// `--sankoff`'s FASTA and every export built from it (`--tnt`, + /// `--phyg`, `--iqtree`), instead of an ordinary, costed 16th alphabet + /// state (the default). For genome-skim/reduced-representation inputs + /// (coverage often < 1x), non-detection is dominated by sampling + /// failure, not true loss — scoring it as a real state risks grouping + /// genomes by shared undersampling rather than shared ancestry. `?` + /// (not `-`) because `-` still carries gap/indel semantics in these + /// tools; a non-detected family is not an observed deletion. See + /// `docmd/theory/evolutionary_distances.md`, "Locus dropout under + /// incomplete coverage". + #[arg(long)] + pub free_loss: bool, + + /// Exclude genome pairs whose raw SNP ratio exceeds this value from the + /// `p_hat` calibration pooled by `--sankoff` — a pair this close to + /// saturation carries no information about `p_hat` and would bias it + /// upward if pooled in (unlike a low eligible-loci count, which barely + /// moves the pooled estimate either way — see design doc). + #[arg(long, default_value = "0.5")] + pub sankoff_ratio_ceiling: f64, + + /// Also write _sankoff.tnt, a ready-to-run TNT script (`proc + /// ;`) for the same matrix/alignment `--sankoff` computes — + /// recoded to TNT's default xread alphabet (0-9A-F only; TNT rejects + /// the wider IUPAC set `--sankoff`'s own output uses unless `nstates + /// dna` is set, which imposes TNT's own incompatible DNA encoding + /// instead) with integer-scaled costs (TNT's smatrix/cost commands + /// reject decimals). Implies `--sankoff`. + #[arg(long)] + pub tnt: bool, + + /// Also write _sankoff.tcm and _sankoff.pg, a + /// custom-alphabet cost matrix and a ready-to-run PhyG script (`read`/ + /// `search`/`report`) for the same matrix/alignment `--sankoff` + /// computes. Reuses `--sankoff`'s own `_sankoff.fasta` directly — PhyG's + /// `tcm:` alphabet is read from the matrix file itself, so the IUPAC+`0` + /// alphabet needs no recoding here, unlike `--tnt`. Implies `--sankoff`. + #[arg(long)] + pub phyg: bool, + + /// Also write _iqtree.model and _iqtree.fasta, a + /// custom-model file and a matching + /// recoded alignment for genuine maximum-likelihood inference with + /// IQ-TREE (`iqtree3 -s ... --seqtype MORPH -m ...+ASC`) — real branch + /// lengths, unlike `--tnt`/`--phyg`'s parsimony step counts. The model + /// is the reversible `Q(i,j) = R(i,j)·π_j` construction: `R` + /// (exchangeability, symmetric) recovered from the same calibrated + /// cost matrix `--sankoff` computes, `π` the real empirical state + /// frequencies counted from the alignment (not IQ-TREE's `+FO`/`+F` — + /// neither applies to a custom-file model, see + /// `docmd/theory/evolutionary_distances.md`). Only the states that + /// actually occur in this alignment are kept, compactly renumbered + /// (IQ-TREE infers its state count from the alignment itself, and a + /// gap in the numbering would silently misalign the model file). + /// Implies `--sankoff`. + #[arg(long)] + pub iqtree: bool, + + /// Scale factor applied before rounding real-valued costs to the + /// integers both `--tnt`'s smatrix/cost commands and `--phyg`'s `tcm:` + /// matrix require. Keep this small: the total tree score is this scale + /// times the sum of per-character costs across every character (908k+ + /// for a typical run here), and there are hints in TNT's own manual + /// that at least some of its internal accumulators are 32-bit — a large + /// scale risks a silent integer overflow (undetectable, not just a + /// crash) far more costly than the resolution a bigger factor would + /// buy. Shared between `--tnt` and `--phyg` rather than split into two + /// flags: both scale the same calibrated matrix for the same reason + /// (integer-only cost commands), and no PhyG-specific accumulator-width + /// constraint has actually been found to justify a different default. + #[arg(long, default_value = "100")] + pub sankoff_cost_scale: f64, + + /// Output prefix: _dist.csv, _shared.csv, + /// _siblings.csv, _rawsnp.csv, _rawsnp_counts.csv, + /// _snp.fasta, _family_overlap.csv, + /// _sankoff_matrix.csv, _sankoff_params.yaml, + /// _sankoff.fasta, _sankoff.tnt, _sankoff.tcm, + /// _sankoff.pg, _iqtree.model, _iqtree.fasta, + /// _nj.nwk, _upgma.nwk. + /// If omitted, the distance matrix is written to stdout. + #[arg(short, long)] + pub output: Option, +} diff --git a/src/obikmer/src/cmd/phylo/family_overlap.rs b/src/obikmer/src/cmd/phylo/family_overlap.rs new file mode 100644 index 00000000..fa01c8e2 --- /dev/null +++ b/src/obikmer/src/cmd/phylo/family_overlap.rs @@ -0,0 +1,97 @@ +use std::io::{BufWriter, Write}; +use std::path::PathBuf; + +use obikindex::{KmerIndex, SnpAlignment}; +use tracing::info; + +// ── Family overlap: shared-family counts and the `--min-shared-family` / +// `--family-overlap` diagnostics built from them ──────────────────────────── +// +// Same variable-family columns as `--snp`'s pseudo-alignment. Off-diagonal +// `[i][j]`: number of columns where both genome `i` and genome `j` carry a +// call (neither is `∅`) — how much informative family content two genomes +// actually share, the direct diagnostic for the rogue-taxon placement seen +// under `--free-loss` (a genome with little overlap with anything else has +// almost nothing left to constrain it). Diagonal `[i][i]` kept, deliberately +// not skipped: with `i == j` the condition "both non-`∅`" degenerates to +// "genome `i` non-`∅`", i.e. the total number of variable families genome +// `i` carries at all — a genome-level count worth having alongside the +// pairwise ones, not a separate computation. + +/// `counts[i][j]` = number of variable-family columns where both genome `i` +/// and genome `j` carry a call (neither is `∅`). Shared between +/// `write_family_overlap_csv` and `--min-shared-family`'s auto-exclusion so +/// both read off the same definition of "shared family". +fn family_overlap_counts(alignment: &SnpAlignment) -> Vec> { + let n = alignment.sequences.len(); + let mut counts = vec![vec![0u64; n]; n]; + for i in 0..n { + for j in 0..n { + counts[i][j] = alignment.sequences[i].iter().zip(alignment.sequences[j].iter()) + .filter(|&(&a, &b)| a != b'-' && b != b'-') + .count() as u64; + } + } + counts +} + +/// Mean of row `i` in a `family_overlap_counts` matrix, excluding the +/// diagonal — how much informative content genome `i` shares with the +/// *average* other genome, the statistic `--min-shared-family` thresholds. +fn mean_offdiag(counts: &[Vec], i: usize) -> f64 { + let n = counts.len(); + let sum: u64 = (0..n).filter(|&j| j != i).map(|j| counts[i][j]).sum(); + sum as f64 / (n - 1) as f64 +} + +pub(super) fn write_family_overlap_csv(alignment: &SnpAlignment, labels: &[String], output: &Option) { + let path = output.as_ref() + .map(|p| format!("{}_family_overlap.csv", p.display())) + .unwrap_or_else(|| "family_overlap.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); + })); + let n = labels.len(); + let counts = family_overlap_counts(alignment); + write!(f, "genome").unwrap(); + for g in labels { write!(f, ",{g}").unwrap(); } + writeln!(f).unwrap(); + for (i, gi) in labels.iter().enumerate() { + write!(f, "{gi}").unwrap(); + for j in 0..n { + write!(f, ",{}", counts[i][j]).unwrap(); + } + writeln!(f).unwrap(); + } + info!("family overlap matrix → {path}"); +} + +/// Sets `mask[i] = true` for every genome whose mean shared-family count +/// (`mean_offdiag`) falls below `threshold`, skipping genomes already +/// excluded (`mask[i]` already `true`, e.g. via `--exclude-genome`). Builds +/// its own `SnpAlignment` pass — same redundant-per-flag pattern already +/// used throughout `run()` (`--snp`/`--sankoff`/`--family-overlap` each call +/// `snp_pseudo_alignment` independently too). +pub(super) fn apply_min_shared_family_exclusion( + idx: &KmerIndex, + labels: &[String], + threshold: f64, + mask: &mut [bool], +) { + let alignment = idx.snp_pseudo_alignment().unwrap_or_else(|e| { + eprintln!("error computing SNP pseudo-alignment for --min-shared-family: {e}"); + std::process::exit(1); + }); + let counts = family_overlap_counts(&alignment); + for (i, label) in labels.iter().enumerate() { + if mask[i] { + continue; // already excluded via --exclude-genome + } + let mean = mean_offdiag(&counts, i); + if mean < threshold { + info!("--min-shared-family: excluding {label} (mean shared families = {mean:.1} < {threshold})"); + mask[i] = true; + } + } +} diff --git a/src/obikmer/src/cmd/phylo/iqtree.rs b/src/obikmer/src/cmd/phylo/iqtree.rs index 0a1ac655..cd906841 100644 --- a/src/obikmer/src/cmd/phylo/iqtree.rs +++ b/src/obikmer/src/cmd/phylo/iqtree.rs @@ -68,13 +68,59 @@ impl CompactAlphabet { } } -fn compact_alphabet(alignment: &SnpAlignment) -> CompactAlphabet { +/// Under `--free-loss`, non-detection (`-`) becomes IQ-TREE's own missing +/// symbol (`?`) — ignored when IQ-TREE checks a site's constancy for +/// `+ASC`. A family kept as "variable" by `snp_pseudo_alignment` +/// (`family_size() >= 2`, a whole-annex property, oblivious to any one +/// column's actual calls) can still turn constant *among the genomes that +/// actually have data* once the non-detected ones are excluded from that +/// check — the same failure mode as the `--exclude-genome`/`drop_excluded` +/// fix in `mod.rs` (see `docmd/theory/evolutionary_distances.md`, "Two +/// consistency bugs found and fixed post-implementation"), just triggered +/// by hiding cells instead of dropping whole rows. Same remedy: rescan +/// columns treating `-` as ignored, drop any where the remaining calls +/// agree on a single state. Parsimony (`--tnt`/`--phyg`) has no +/// no-invariant-site requirement, so this only runs on IQ-TREE's own copy +/// of the alignment, never mutating the one the caller also hands to those +/// two exports. +fn drop_ascertainment_noninformative(alignment: &SnpAlignment) -> SnpAlignment { + let n_sites = alignment.sequences.first().map(|s| s.len()).unwrap_or(0); + let keep: Vec = (0..n_sites).map(|site| { + let mut first: Option = None; + for seq in &alignment.sequences { + let b = seq[site]; + if b == b'-' { + continue; + } + match first { + None => first = Some(b), + Some(f) if f != b => return true, + _ => {} + } + } + false // all calls missing, or all calls agree — non-informative + }).collect(); + + let sequences = alignment.sequences.iter() + .map(|seq| seq.iter().zip(keep.iter()).filter(|&(_, &k)| k).map(|(&b, _)| b).collect()) + .collect(); + SnpAlignment { sequences } +} + +fn compact_alphabet(alignment: &SnpAlignment, free_loss: bool) -> CompactAlphabet { let iupac_to_state = state_index_table(); let mut occurs = [false; 16]; let mut counts = [0u64; 16]; for seq in &alignment.sequences { for &b in seq { + if free_loss && b == b'-' { + // `?`: IQ-TREE's own missing-data symbol for `--seqtype + // MORPH`, marginalised by Felsenstein pruning — not a + // numbered state, so excluded from `occurs`/`counts` and + // from the compact alphabet built below. + continue; + } let b = if b == b'-' { b'0' } else { b }; let state = iupac_to_state[b as usize] as usize; occurs[state] = true; @@ -133,6 +179,7 @@ fn write_iqtree_alignment( labels: &[String], alphabet: &CompactAlphabet, output: &Option, + free_loss: bool, ) -> (String, usize) { let iupac_to_state = state_index_table(); @@ -146,6 +193,9 @@ fn write_iqtree_alignment( 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 free_loss && b == b'-' { + return b'?'; + } let b = if b == b'-' { b'0' } else { b }; let old = iupac_to_state[b as usize] as usize; let compact = alphabet.old_to_compact[old] @@ -165,10 +215,27 @@ pub(super) fn write_iqtree( alignment: &SnpAlignment, labels: &[String], output: &Option, + free_loss: bool, ) { - let alphabet = compact_alphabet(alignment); + let filtered; + let alignment = if free_loss { + let before = alignment.sequences.first().map(|s| s.len()).unwrap_or(0); + filtered = drop_ascertainment_noninformative(alignment); + let after = filtered.sequences.first().map(|s| s.len()).unwrap_or(0); + if after != before { + info!( + "--free-loss: {before} → {after} sites (dropped columns non-informative once `-` \ + is treated as missing — required for +ASC)" + ); + } + &filtered + } else { + alignment + }; + + let alphabet = compact_alphabet(alignment, free_loss); let model_path = write_iqtree_model(matrix, &alphabet, output); - let (fasta_path, n_sites) = write_iqtree_alignment(alignment, labels, &alphabet, output); + let (fasta_path, n_sites) = write_iqtree_alignment(alignment, labels, &alphabet, output, free_loss); let prefix_name = output.as_ref() .and_then(|p| p.file_name()) diff --git a/src/obikmer/src/cmd/phylo/mod.rs b/src/obikmer/src/cmd/phylo/mod.rs index a25f7a39..1f08cdc5 100644 --- a/src/obikmer/src/cmd/phylo/mod.rs +++ b/src/obikmer/src/cmd/phylo/mod.rs @@ -1,225 +1,29 @@ +mod args; +mod family_overlap; mod iqtree; +mod outputs; mod phyg; mod sankoff; mod tnt; use std::io::{self, BufWriter, Write}; -use std::path::PathBuf; -use clap::Args; use kodama::{Method, linkage}; -use obifastwrite::{JsonVal, write_record}; use obikindex::{ - DistanceMetric, KmerIndex, RawSnpDistanceOutput, - SiblingAnnexStats, SnpAlignment, + KmerIndex, RawSnpDistanceOutput, SnpAlignment, cardinality_transition_probs, composition_transition_probs, pairwise_cost_matrix, }; use speedytree::{DistanceMatrix, Hybrid, NeighborJoiningSolver, to_newick}; use tracing::info; +pub use args::PhyloArgs; +use family_overlap::{apply_min_shared_family_exclusion, write_family_overlap_csv}; use iqtree::write_iqtree; +use outputs::{write_raw_snp_counts_csv, write_raw_snp_distance_csv, write_sibling_stats_csv, write_snp_fasta, upgma_to_newick}; 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, - Mash, - Hamming, - BrayCurtis, - #[value(name = "relfreq-bray-curtis")] - RelfreqBrayCurtis, - Euclidean, - #[value(name = "relfreq-euclidean")] - RelfreqEuclidean, - Hellinger, - #[value(name = "hellinger-euclidean")] - HellingerEuclidean, -} - -impl From for DistanceMetric { - fn from(m: MetricArg) -> Self { - match m { - MetricArg::Jaccard => DistanceMetric::Jaccard, - MetricArg::Mash => DistanceMetric::Mash, - MetricArg::Hamming => DistanceMetric::Hamming, - MetricArg::BrayCurtis => DistanceMetric::BrayCurtis, - MetricArg::RelfreqBrayCurtis => DistanceMetric::RelfreqBrayCurtis, - MetricArg::Euclidean => DistanceMetric::Euclidean, - MetricArg::RelfreqEuclidean => DistanceMetric::RelfreqEuclidean, - MetricArg::Hellinger => DistanceMetric::Hellinger, - MetricArg::HellingerEuclidean => DistanceMetric::HellingerEuclidean, - } - } -} - -#[derive(Args)] -pub struct PhyloArgs { - /// Index directory - pub index: PathBuf, - - /// Distance metric to compute - #[arg(long, value_enum, default_value = "jaccard")] - pub metric: MetricArg, - - /// Minimum count to consider a kmer present when computing Jaccard on count indexes - #[arg(long, default_value = "1")] - pub presence_threshold: u32, - - /// Also output the shared-kmer count matrix (CSV) - #[arg(long)] - pub shared_kmers: bool, - - /// Compute and write a Neighbor-Joining tree (Newick) - #[arg(long)] - pub nj: bool, - - /// Compute and write a UPGMA tree (Newick) - #[arg(long)] - pub upgma: bool, - - /// Build the sibling-count/minorant annex on this (multi-genome) index - /// — see `docmd/theory/evolutionary_distances.md`, Step 2b. Construction - /// only; does not by itself compute or write any statistics. - #[arg(long)] - pub sibling_annex: bool, - - /// Exclude a genome (by its exact label) from every computation below - /// that reads the sibling annex — `--raw-snp-distance`/`--raw-snp-counts`, - /// `--snp`, and `--sankoff` (and everything `--sankoff` implies: the - /// cardinality/composition transition models, the exported - /// matrix/alignment, `--tnt`/`--phyg`/`--iqtree`). Repeatable. Does - /// *not* affect the plain `--metric` distance matrix/NJ/UPGMA path (a - /// different, unrelated computation). Applied by zeroing the excluded - /// genome's row/column after `raw_snp_distance` runs (a pair with zero - /// counts is already skipped by `base_pair_tally`/`cardinality_tally`, - /// so this needs no change to the underlying traversal) and by - /// dropping its row from `snp_pseudo_alignment`'s output — the annex - /// is still built/scanned for the excluded genome too, just not used - /// afterward. For a genome with almost no informative sites shared - /// with anything else (see `docmd/theory/evolutionary_distances.md`, - /// the IQ-TREE/Mash rogue-taxon discussion), its presence can - /// otherwise silently bias the transition models. - #[arg(long = "exclude-genome", value_name = "LABEL")] - pub exclude_genome: Vec, - - /// Tally the sibling-count distribution (CSV) of an already-built annex - /// (run with `--sibling-annex` first, in this invocation or an earlier - /// one). A separate, occasional diagnostic pass — not run every time the - /// annex itself is (re)built. - #[arg(long)] - pub sibling_stats: bool, - - /// Compute the raw p-distance restricted to loci that are single-copy - /// in both genomes of each pair (an already-built sibling annex is - /// required — run with `--sibling-annex` first, in this invocation or - /// an earlier one). A quick way to test the central-position SNP - /// estimator against a real index; not the full `SnpTally` design. - #[arg(long)] - pub raw_snp_distance: bool, - - /// Write the raw per-pair counts (`n_snp`, `n_shared`, `n_eligible`) - /// behind `--raw-snp-distance`'s ratio, one row per genome pair — a - /// diagnostic table, not a matrix. The ratio alone can't distinguish - /// "identical at every eligible locus" from "almost no eligible loci - /// at all" (e.g. `0.0` from 0/2 looks the same as `0.0` from 0/2000), - /// and that distinction matters a lot for genome pairs near the edge - /// of what central-position families can resolve (see - /// `docmd/theory/evolutionary_distances.md`, "Run 3" and the - /// IQ-TREE/Mash comparison). Same annex requirement as - /// `--raw-snp-distance`. - #[arg(long)] - pub raw_snp_counts: bool, - - /// Write a SNP-only pseudo-alignment (FASTA, IUPAC-coded) from an - /// already-built sibling annex — one row per genome, one column per - /// variable family (monomorphic families skipped), no flanking - /// sequence. See `docmd/theory/evolutionary_distances.md`, - /// "Multi-genome framing: family as pseudo-alignment column". - #[arg(long)] - pub snp: bool, - - /// Calibrate a 16-state Sankoff cost matrix and its matching - /// pseudo-alignment from an already-built sibling annex (run with - /// `--sibling-annex` first, in this invocation or an earlier one), for - /// use with TNT/PhyG. See `docmd/theory/evolutionary_distances.md`, - /// "Sankoff parsimony as the resolution of the 16-state model problem". - #[arg(long)] - pub sankoff: bool, - - /// Exclude genome pairs whose raw SNP ratio exceeds this value from the - /// `p_hat` calibration pooled by `--sankoff` — a pair this close to - /// saturation carries no information about `p_hat` and would bias it - /// upward if pooled in (unlike a low eligible-loci count, which barely - /// moves the pooled estimate either way — see design doc). - #[arg(long, default_value = "0.5")] - pub sankoff_ratio_ceiling: f64, - - /// Also write _sankoff.tnt, a ready-to-run TNT script (`proc - /// ;`) for the same matrix/alignment `--sankoff` computes — - /// recoded to TNT's default xread alphabet (0-9A-F only; TNT rejects - /// the wider IUPAC set `--sankoff`'s own output uses unless `nstates - /// dna` is set, which imposes TNT's own incompatible DNA encoding - /// instead) with integer-scaled costs (TNT's smatrix/cost commands - /// reject decimals). Implies `--sankoff`. - #[arg(long)] - pub tnt: bool, - - /// Also write _sankoff.tcm and _sankoff.pg, a - /// custom-alphabet cost matrix and a ready-to-run PhyG script (`read`/ - /// `search`/`report`) for the same matrix/alignment `--sankoff` - /// computes. Reuses `--sankoff`'s own `_sankoff.fasta` directly — PhyG's - /// `tcm:` alphabet is read from the matrix file itself, so the IUPAC+`0` - /// alphabet needs no recoding here, unlike `--tnt`. Implies `--sankoff`. - #[arg(long)] - pub phyg: bool, - - /// Also write _iqtree.model and _iqtree.fasta, a - /// custom-model file and a matching - /// recoded alignment for genuine maximum-likelihood inference with - /// IQ-TREE (`iqtree3 -s ... --seqtype MORPH -m ...+ASC`) — real branch - /// lengths, unlike `--tnt`/`--phyg`'s parsimony step counts. The model - /// is the reversible `Q(i,j) = R(i,j)·π_j` construction: `R` - /// (exchangeability, symmetric) recovered from the same calibrated - /// cost matrix `--sankoff` computes, `π` the real empirical state - /// frequencies counted from the alignment (not IQ-TREE's `+FO`/`+F` — - /// neither applies to a custom-file model, see - /// `docmd/theory/evolutionary_distances.md`). Only the states that - /// actually occur in this alignment are kept, compactly renumbered - /// (IQ-TREE infers its state count from the alignment itself, and a - /// gap in the numbering would silently misalign the model file). - /// Implies `--sankoff`. - #[arg(long)] - pub iqtree: bool, - - /// Scale factor applied before rounding real-valued costs to the - /// integers both `--tnt`'s smatrix/cost commands and `--phyg`'s `tcm:` - /// matrix require. Keep this small: the total tree score is this scale - /// times the sum of per-character costs across every character (908k+ - /// for a typical run here), and there are hints in TNT's own manual - /// that at least some of its internal accumulators are 32-bit — a large - /// scale risks a silent integer overflow (undetectable, not just a - /// crash) far more costly than the resolution a bigger factor would - /// buy. Shared between `--tnt` and `--phyg` rather than split into two - /// flags: both scale the same calibrated matrix for the same reason - /// (integer-only cost commands), and no PhyG-specific accumulator-width - /// constraint has actually been found to justify a different default. - #[arg(long, default_value = "100")] - pub sankoff_cost_scale: f64, - - /// Output prefix: _dist.csv, _shared.csv, - /// _siblings.csv, _rawsnp.csv, _rawsnp_counts.csv, - /// _snp.fasta, - /// _sankoff_matrix.csv, _sankoff_params.yaml, - /// _sankoff.fasta, _sankoff.tnt, _sankoff.tcm, - /// _sankoff.pg, _iqtree.model, _iqtree.fasta, - /// _nj.nwk, _upgma.nwk. - /// If omitted, the distance matrix is written to stdout. - #[arg(short, long)] - pub output: Option, -} - pub fn run(args: PhyloArgs) { let idx = KmerIndex::open(&args.index).unwrap_or_else(|e| { eprintln!("error opening index: {e}"); @@ -250,6 +54,9 @@ pub fn run(args: PhyloArgs) { } } } + if let Some(threshold) = args.min_shared_family { + apply_min_shared_family_exclusion(&idx, &labels, threshold, &mut mask); + } mask }; let zero_excluded_pairs = |result: &mut RawSnpDistanceOutput| { @@ -352,6 +159,14 @@ pub fn run(args: PhyloArgs) { let (alignment, kept_labels) = drop_excluded(alignment); write_snp_fasta(&alignment, &kept_labels, &args.output); } + if args.family_overlap { + let alignment = idx.snp_pseudo_alignment().unwrap_or_else(|e| { + eprintln!("error computing SNP pseudo-alignment: {e}"); + std::process::exit(1); + }); + let (alignment, kept_labels) = drop_excluded(alignment); + write_family_overlap_csv(&alignment, &kept_labels, &args.output); + } if args.sankoff || args.tnt || args.phyg || args.iqtree { let mut raw = idx.raw_snp_distance().unwrap_or_else(|e| { eprintln!("error computing raw SNP distance: {e}"); @@ -378,16 +193,16 @@ pub fn run(args: PhyloArgs) { std::process::exit(1); }); let (alignment, kept_labels) = drop_excluded(alignment); - write_sankoff_alignment_fasta(&alignment, &kept_labels, &args.output); + write_sankoff_alignment_fasta(&alignment, &kept_labels, &args.output, args.free_loss); if args.tnt { - write_sankoff_tnt(&matrix, &alignment, &kept_labels, &args.output, args.sankoff_cost_scale); + write_sankoff_tnt(&matrix, &alignment, &kept_labels, &args.output, args.sankoff_cost_scale, args.free_loss); } if args.phyg { write_sankoff_phyg(&matrix, &args.output, args.sankoff_cost_scale); } if args.iqtree { - write_iqtree(&matrix, &alignment, &kept_labels, &args.output); + write_iqtree(&matrix, &alignment, &kept_labels, &args.output, args.free_loss); } } @@ -403,6 +218,7 @@ pub fn run(args: PhyloArgs) { || args.raw_snp_distance || args.raw_snp_counts || args.snp + || args.family_overlap || args.sankoff || args.tnt || args.phyg @@ -521,165 +337,3 @@ pub fn run(args: PhyloArgs) { info!("UPGMA tree → {path}"); } } - -// ── Family-size distribution → CSV ────────────────────────────────────────── -// -// Each row is a family (the up-to-4 k-mers sharing flanks, differing only at -// the centre), counted once — at its minorant — regardless of how many of -// its members are observed. Family size 1..4 (not "sibling count" 0..3): -// see `docmd/theory/evolutionary_distances.md`, "Definitions". - -fn write_sibling_stats_csv(stats: &SiblingAnnexStats, labels: &[String], output: &Option) { - // One row per genome (4 columns, family size 1-4: number of families of - // that size for which the genome carries at least one member), plus a - // `global` row — the actual deduplicated family-size histogram - // (`stats.counts`), NOT a sum of the per-genome columns (a family shared - // by several genomes would otherwise be counted once per genome it - // appears in, inflating the total beyond the real family count). - let path = output.as_ref() - .map(|p| format!("{}_siblings.csv", p.display())) - .unwrap_or_else(|| "siblings.csv".into()); - let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| { - eprintln!("error creating {path}: {e}"); - std::process::exit(1); - })); - writeln!(f, "genome,1,2,3,4").unwrap(); - for (label, counts) in labels.iter().zip(stats.per_genome.iter()) { - writeln!(f, "{label},{},{},{},{}", counts[0], counts[1], counts[2], counts[3]).unwrap(); - } - writeln!( - f, "global,{},{},{},{}", - stats.counts[0], stats.counts[1], stats.counts[2], stats.counts[3], - ).unwrap(); - let total: u64 = stats.counts.iter().sum(); - info!("family-size distribution → {path} (total {total} famil{})", - if total == 1 { "y" } else { "ies" }); -} - -// ── Raw single-copy SNP distance → CSV ────────────────────────────────────── -// -// p_hat[i,j] = snp[i,j] / (snp[i,j] + shared[i,j]) over loci single-copy in -// both i and j — see `RawSnpDistanceOutput` / `KmerIndex::raw_snp_distance`. -// A single file: the distance matrix, with an eligible-loci count alongside -// each value so a 0/0 pair (no eligible locus at all) is distinguishable -// from a genuinely identical pair. - -fn write_raw_snp_distance_csv(result: &RawSnpDistanceOutput, labels: &[String], output: &Option) { - let path = output.as_ref() - .map(|p| format!("{}_rawsnp.csv", p.display())) - .unwrap_or_else(|| "rawsnp.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); - })); - let n = labels.len(); - write!(f, "genome").unwrap(); - for g in labels { write!(f, ",{g}").unwrap(); } - writeln!(f).unwrap(); - for (i, g) in labels.iter().enumerate() { - write!(f, "{g}").unwrap(); - for j in 0..n { - let snp = result.snp[[i, j]]; - let shared = result.shared[[i, j]]; - let eligible = snp + shared; - if eligible == 0 { - write!(f, ",NA").unwrap(); - } else { - write!(f, ",{:.6}", snp as f64 / eligible as f64).unwrap(); - } - } - writeln!(f).unwrap(); - } - info!("raw single-copy SNP distance matrix → {path}"); -} - -// ── Raw single-copy SNP distance → per-pair diagnostic counts ────────────── -// -// A pair table (one row per unordered genome pair), not a matrix: the ratio -// alone can't distinguish "identical across every eligible locus" from -// "almost no eligible locus at all" — both can read `0.0`/`NA` in -// `--raw-snp-distance`'s output. Distinguishing them matters most exactly -// where it's easy to miss: genome pairs near the edge of what -// central-position families can resolve at all (deep cross-lineage splits, -// see `docmd/theory/evolutionary_distances.md`, "Run 3" and the later -// IQ-TREE/Mash comparison — a `ratio=0.0` backed by 2 eligible loci is not -// the same claim as one backed by 2000). - -fn write_raw_snp_counts_csv(result: &RawSnpDistanceOutput, labels: &[String], output: &Option) { - let path = output.as_ref() - .map(|p| format!("{}_rawsnp_counts.csv", p.display())) - .unwrap_or_else(|| "rawsnp_counts.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); - })); - let n = labels.len(); - writeln!(f, "genome_a,genome_b,n_snp,n_shared,n_eligible,ratio").unwrap(); - for i in 0..n { - for j in (i + 1)..n { - let snp = result.snp[[i, j]]; - let shared = result.shared[[i, j]]; - let eligible = snp + shared; - write!(f, "{},{},{snp},{shared},{eligible}", labels[i], labels[j]).unwrap(); - if eligible == 0 { - writeln!(f, ",NA").unwrap(); - } else { - writeln!(f, ",{:.6}", snp as f64 / eligible as f64).unwrap(); - } - } - } - info!("raw single-copy SNP distance counts (diagnostic) → {path}"); -} - -// ── SNP-only pseudo-alignment → FASTA ─────────────────────────────────────── -// -// One record per genome, IUPAC-coded, no flanking sequence — see -// `SnpAlignment` / `KmerIndex::snp_pseudo_alignment`. Uses the project's -// existing FASTA writer (`obifastwrite::write_record`) rather than -// hand-rolling one. - -fn write_snp_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option) { - let path = output.as_ref() - .map(|p| format!("{}_snp.fasta", p.display())) - .unwrap_or_else(|| "snp.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()) { - write_record(seq, 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!("SNP pseudo-alignment → {path} ({n_sites} site{})", - if n_sites == 1 { "" } else { "s" }); -} - -// ── UPGMA Newick from kodama dendrogram ─────────────────────────────────────── - -fn upgma_to_newick(dendro: &kodama::Dendrogram, names: &[String]) -> String { - let n = names.len(); - // node_labels[i]: Newick subtree string for node i (leaves 0..n, internals n..) - let mut labels: Vec = names.to_vec(); - // height of each node: leaves = 0, internal = dissimilarity/2 - let mut heights: Vec = vec![0.0; 2 * n - 1]; - - for (k, step) in dendro.steps().iter().enumerate() { - let new_node = n + k; - let h = step.dissimilarity / 2.0; - heights[new_node] = h; - let c1 = step.cluster1; - let c2 = step.cluster2; - let bl1 = (h - heights[c1]).max(0.0); - let bl2 = (h - heights[c2]).max(0.0); - labels.push(format!( - "({label1}:{bl1:.6},{label2}:{bl2:.6})", - label1 = labels[c1], - label2 = labels[c2], - )); - } - - format!("{};", labels.last().unwrap()) -} diff --git a/src/obikmer/src/cmd/phylo/outputs.rs b/src/obikmer/src/cmd/phylo/outputs.rs new file mode 100644 index 00000000..97b548fd --- /dev/null +++ b/src/obikmer/src/cmd/phylo/outputs.rs @@ -0,0 +1,168 @@ +use std::io::{BufWriter, Write}; +use std::path::PathBuf; + +use obifastwrite::{JsonVal, write_record}; +use obikindex::{RawSnpDistanceOutput, SiblingAnnexStats, SnpAlignment}; +use tracing::info; + +// ── Family-size distribution → CSV ────────────────────────────────────────── +// +// Each row is a family (the up-to-4 k-mers sharing flanks, differing only at +// the centre), counted once — at its minorant — regardless of how many of +// its members are observed. Family size 1..4 (not "sibling count" 0..3): +// see `docmd/theory/evolutionary_distances.md`, "Definitions". + +pub(super) fn write_sibling_stats_csv(stats: &SiblingAnnexStats, labels: &[String], output: &Option) { + // One row per genome (4 columns, family size 1-4: number of families of + // that size for which the genome carries at least one member), plus a + // `global` row — the actual deduplicated family-size histogram + // (`stats.counts`), NOT a sum of the per-genome columns (a family shared + // by several genomes would otherwise be counted once per genome it + // appears in, inflating the total beyond the real family count). + let path = output.as_ref() + .map(|p| format!("{}_siblings.csv", p.display())) + .unwrap_or_else(|| "siblings.csv".into()); + let mut f = BufWriter::new(std::fs::File::create(&path).unwrap_or_else(|e| { + eprintln!("error creating {path}: {e}"); + std::process::exit(1); + })); + writeln!(f, "genome,1,2,3,4").unwrap(); + for (label, counts) in labels.iter().zip(stats.per_genome.iter()) { + writeln!(f, "{label},{},{},{},{}", counts[0], counts[1], counts[2], counts[3]).unwrap(); + } + writeln!( + f, "global,{},{},{},{}", + stats.counts[0], stats.counts[1], stats.counts[2], stats.counts[3], + ).unwrap(); + let total: u64 = stats.counts.iter().sum(); + info!("family-size distribution → {path} (total {total} famil{})", + if total == 1 { "y" } else { "ies" }); +} + +// ── Raw single-copy SNP distance → CSV ────────────────────────────────────── +// +// p_hat[i,j] = snp[i,j] / (snp[i,j] + shared[i,j]) over loci single-copy in +// both i and j — see `RawSnpDistanceOutput` / `KmerIndex::raw_snp_distance`. +// A single file: the distance matrix, with an eligible-loci count alongside +// each value so a 0/0 pair (no eligible locus at all) is distinguishable +// from a genuinely identical pair. + +pub(super) fn write_raw_snp_distance_csv(result: &RawSnpDistanceOutput, labels: &[String], output: &Option) { + let path = output.as_ref() + .map(|p| format!("{}_rawsnp.csv", p.display())) + .unwrap_or_else(|| "rawsnp.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); + })); + let n = labels.len(); + write!(f, "genome").unwrap(); + for g in labels { write!(f, ",{g}").unwrap(); } + writeln!(f).unwrap(); + for (i, g) in labels.iter().enumerate() { + write!(f, "{g}").unwrap(); + for j in 0..n { + let snp = result.snp[[i, j]]; + let shared = result.shared[[i, j]]; + let eligible = snp + shared; + if eligible == 0 { + write!(f, ",NA").unwrap(); + } else { + write!(f, ",{:.6}", snp as f64 / eligible as f64).unwrap(); + } + } + writeln!(f).unwrap(); + } + info!("raw single-copy SNP distance matrix → {path}"); +} + +// ── Raw single-copy SNP distance → per-pair diagnostic counts ────────────── +// +// A pair table (one row per unordered genome pair), not a matrix: the ratio +// alone can't distinguish "identical across every eligible locus" from +// "almost no eligible locus at all" — both can read `0.0`/`NA` in +// `--raw-snp-distance`'s output. Distinguishing them matters most exactly +// where it's easy to miss: genome pairs near the edge of what +// central-position families can resolve at all (deep cross-lineage splits, +// see `docmd/theory/evolutionary_distances.md`, "Run 3" and the later +// IQ-TREE/Mash comparison — a `ratio=0.0` backed by 2 eligible loci is not +// the same claim as one backed by 2000). + +pub(super) fn write_raw_snp_counts_csv(result: &RawSnpDistanceOutput, labels: &[String], output: &Option) { + let path = output.as_ref() + .map(|p| format!("{}_rawsnp_counts.csv", p.display())) + .unwrap_or_else(|| "rawsnp_counts.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); + })); + let n = labels.len(); + writeln!(f, "genome_a,genome_b,n_snp,n_shared,n_eligible,ratio").unwrap(); + for i in 0..n { + for j in (i + 1)..n { + let snp = result.snp[[i, j]]; + let shared = result.shared[[i, j]]; + let eligible = snp + shared; + write!(f, "{},{},{snp},{shared},{eligible}", labels[i], labels[j]).unwrap(); + if eligible == 0 { + writeln!(f, ",NA").unwrap(); + } else { + writeln!(f, ",{:.6}", snp as f64 / eligible as f64).unwrap(); + } + } + } + info!("raw single-copy SNP distance counts (diagnostic) → {path}"); +} + +// ── SNP-only pseudo-alignment → FASTA ─────────────────────────────────────── +// +// One record per genome, IUPAC-coded, no flanking sequence — see +// `SnpAlignment` / `KmerIndex::snp_pseudo_alignment`. Uses the project's +// existing FASTA writer (`obifastwrite::write_record`) rather than +// hand-rolling one. + +pub(super) fn write_snp_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option) { + let path = output.as_ref() + .map(|p| format!("{}_snp.fasta", p.display())) + .unwrap_or_else(|| "snp.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()) { + write_record(seq, 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!("SNP pseudo-alignment → {path} ({n_sites} site{})", + if n_sites == 1 { "" } else { "s" }); +} + +// ── UPGMA Newick from kodama dendrogram ─────────────────────────────────────── + +pub(super) fn upgma_to_newick(dendro: &kodama::Dendrogram, names: &[String]) -> String { + let n = names.len(); + // node_labels[i]: Newick subtree string for node i (leaves 0..n, internals n..) + let mut labels: Vec = names.to_vec(); + // height of each node: leaves = 0, internal = dissimilarity/2 + let mut heights: Vec = vec![0.0; 2 * n - 1]; + + for (k, step) in dendro.steps().iter().enumerate() { + let new_node = n + k; + let h = step.dissimilarity / 2.0; + heights[new_node] = h; + let c1 = step.cluster1; + let c2 = step.cluster2; + let bl1 = (h - heights[c1]).max(0.0); + let bl2 = (h - heights[c2]).max(0.0); + labels.push(format!( + "({label1}:{bl1:.6},{label2}:{bl2:.6})", + label1 = labels[c1], + label2 = labels[c2], + )); + } + + format!("{};", labels.last().unwrap()) +} diff --git a/src/obikmer/src/cmd/phylo/sankoff.rs b/src/obikmer/src/cmd/phylo/sankoff.rs index c8faaf2b..ab96e7d3 100644 --- a/src/obikmer/src/cmd/phylo/sankoff.rs +++ b/src/obikmer/src/cmd/phylo/sankoff.rs @@ -11,9 +11,15 @@ use tracing::info; // `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. +// rather than our "family absent" state. Unless `free_loss` (`--free-loss`) +// is set, in which case `∅` is recoded to `?` instead — TNT/PhyG's own +// missing-data symbol, deliberately *not* `-` (still gap/indel semantics in +// both tools) — so non-detection costs nothing rather than being scored as +// an ordinary, calibrated state transition. See +// `docmd/theory/evolutionary_distances.md`, "Locus dropout under incomplete +// coverage". -pub(super) fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option) { +pub(super) fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option, free_loss: bool) { let path = output.as_ref() .map(|p| format!("{}_sankoff.fasta", p.display())) .unwrap_or_else(|| "sankoff.fasta".into()); @@ -21,9 +27,10 @@ pub(super) fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[ eprintln!("error creating {path}: {e}"); std::process::exit(1); })); + let absent_symbol = if free_loss { b'?' } else { b'0' }; 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(); + let recoded: Vec = seq.iter().map(|&b| if b == b'-' { absent_symbol } 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); diff --git a/src/obikmer/src/cmd/phylo/tnt.rs b/src/obikmer/src/cmd/phylo/tnt.rs index 48f5ad53..aa16b419 100644 --- a/src/obikmer/src/cmd/phylo/tnt.rs +++ b/src/obikmer/src/cmd/phylo/tnt.rs @@ -28,6 +28,7 @@ pub(super) fn write_sankoff_tnt( labels: &[String], output: &Option, cost_scale: f64, + free_loss: bool, ) { let path = output.as_ref() .map(|p| format!("{}_sankoff.tnt", p.display())) @@ -49,6 +50,13 @@ pub(super) fn write_sankoff_tnt( for (label, seq) in labels.iter().zip(alignment.sequences.iter()) { write!(f, "{label} ").unwrap(); for &b in seq { + if free_loss && b == b'-' { + // `?`: TNT's own missing-data symbol, read directly, not + // routed through `TNT_STATE_SYMBOL` (there is no state for + // it) — see `write_sankoff_alignment_fasta`'s doc comment. + write!(f, "?").unwrap(); + continue; + } 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();