Push zpwxxpnpktps #67
@@ -25,6 +25,7 @@ All functionality is exposed through a single binary, `obikmer`, organized as su
|
||||
| [`dump`](usage/dump.md) | Dump indexed kmers as CSV |
|
||||
| [`annotate`](usage/annotate.md) | Add, update, or dump genome metadata |
|
||||
| [`phylo`](usage/phylo.md) | Compute pairwise evolutionary-distance proxies, trees, and phylogenetic exports |
|
||||
| [`name-tree`](usage/name-tree.md) | Translate a TNT/PhyG numeric-label tree export back to real taxon names |
|
||||
| [`unitig`](usage/unitig.md) | Dump the unitigs of an index as FASTA |
|
||||
| [`estimate`](usage/estimate.md) | Estimate approximate-index parameters before indexing |
|
||||
| [`reindex`](usage/reindex.md) | Convert an index's evidence representation (exact ↔ approximate) |
|
||||
|
||||
@@ -0,0 +1,21 @@
|
||||
# name-tree
|
||||
|
||||
Translate a numerically-labelled tree export (TNT, PhyG, or any plain Newick file with bare `1`, `2`, `3`, … leaf labels) back to real taxon names, reading the label order from the FASTA that produced it.
|
||||
|
||||
```bash
|
||||
obikmer name-tree TREE --fasta FASTA -o OUTPUT
|
||||
```
|
||||
|
||||
## Arguments
|
||||
|
||||
| Argument | Description |
|
||||
|---|---|
|
||||
| `TREE` | Tree file to translate — a TNT-style NEXUS export (`tree NAME = [&U] ...;`) or a plain Newick file |
|
||||
| `--fasta` | FASTA file whose record order gives the numeric taxon labels (1-based) — typically the `_sankoff.fasta`/`_snp.fasta` used to produce `TREE` |
|
||||
| `-o, --output` | Output NEXUS file path |
|
||||
|
||||
## Output
|
||||
|
||||
A NEXUS file with a `taxa` block, a `translate` table (numeric label → taxon name, from `--fasta`'s header order), and every tree found in `TREE`, topology unchanged — readable directly in FigTree, PearTree, `ape` (R), etc.
|
||||
|
||||
`--tnt`'s and `--phyg`'s exports (see [phylo](phylo.md)) both number taxa `1..N` in the same order as the pseudo-alignment FASTA they were built from (`<prefix>_sankoff.fasta`), so pass that same file as `--fasta` here.
|
||||
@@ -93,7 +93,7 @@ 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.
|
||||
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. `--free-loss` also drops the cardinality-transition cost between any two states, not just to/from the empty one: whether a genome shows 1 vs. 2 (etc.) detected members of a family it does carry is exactly as vulnerable to sampling failure as whether the family was detected at all, so gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family. 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
|
||||
|
||||
@@ -103,6 +103,8 @@ All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment,
|
||||
- **`--phyg`**: a custom cost-matrix file plus a PhyG script reusing the `--sankoff` alignment directly.
|
||||
- **`--iqtree`**: a custom substitution-model file (exchangeability matrix recovered as $R(a,b) = e^{-\text{cost}(a,b)}$, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths (unlike the parsimony step-counts from TNT/PhyG). Only states actually occurring in the alignment are kept and compactly renumbered.
|
||||
|
||||
TNT and PhyG both write trees with bare numeric leaf labels (`1`, `2`, …, in the same order as `<prefix>_sankoff.fasta`). Use [`name-tree`](name-tree.md) on the tool's own tree output plus that same FASTA to get a NEXUS file with real taxon names.
|
||||
|
||||
## Output files
|
||||
|
||||
With `-o/--output PREFIX`, the relevant subset of the following files is written: `<prefix>_dist.csv`, `<prefix>_shared.csv`, `<prefix>_nj.nwk`, `<prefix>_upgma.nwk`, `<prefix>_siblings.csv`, `<prefix>_rawsnp.csv`, `<prefix>_rawsnp_counts.csv`, `<prefix>_snp.fasta`, `<prefix>_family_overlap.csv`, `<prefix>_sankoff_matrix.csv`, `<prefix>_sankoff_params.yaml`, `<prefix>_sankoff.fasta`, `<prefix>_sankoff.tnt`, `<prefix>_sankoff.tcm`, `<prefix>_sankoff.pg`, `<prefix>_iqtree.model`, `<prefix>_iqtree.fasta`. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout.
|
||||
|
||||
@@ -1956,6 +1956,79 @@ 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.
|
||||
|
||||
### `--family-overlap`/`--min-shared-family`, and a second `free_loss` bug: cardinality-transition cost not zeroed (2026-08-13)
|
||||
|
||||
**Diagnostic added: `--family-overlap`.** Real per-genome-pair count of
|
||||
shared variable families (columns where neither genome is `∅`), diagonal
|
||||
kept (degenerates to "genome's own total family count", not a separate
|
||||
computation — see `obikmer/src/cmd/phylo/family_overlap.rs`). On the
|
||||
20-genome benchmark this predicted the `--tnt`/`--iqtree` rogue-placement
|
||||
pattern exactly: Candidozyma shares **0** families with all 19 others (hence
|
||||
its ~0-length branch, placed anywhere); the mean off-diagonal shared-family
|
||||
count per genome splits cleanly below/above ~1000 (`Candidozyma` 0,
|
||||
`Saccharolobus` 0.9, `Opitutus` 20.5, `Wolbachia` 37.1, `Acidobacterium`
|
||||
41.2, `Shouchella` 49.4, `Bacillus` 54.5, `Proteus` 629, `Yersinia` 907.5 —
|
||||
all below 1000 — vs. `Klebsiella` ~5000-6000, `Salmonella` ~8000-17000,
|
||||
`E. coli` ~19000-24000, all comfortably above). `--min-shared-family N`
|
||||
automates exclusion below that threshold, reusing the `--exclude-genome`
|
||||
machinery. Restricting to the 3 genera above threshold (11 genomes: 4
|
||||
*E. coli*, 3 *Klebsiella*, 4 *Salmonella*) gives 3 cleanly separated
|
||||
monophyletic clades in both TNT and IQ-TREE (without `+ASC`, see below) —
|
||||
confirms the family-overlap statistic, not just the tree topology, as the
|
||||
right lens for this failure mode. `N=1000` is benchmark-specific (huge
|
||||
cross-domain divergence), not a universal constant — recompute
|
||||
`--family-overlap` and read its own gap on any other dataset.
|
||||
|
||||
**Second `free_loss` bug, found by re-deriving the design rather than by a
|
||||
failing run: the cardinality-transition cost wasn't zeroed.**
|
||||
`--free-loss` neutralises `∅` (whole-family non-detection) by recoding it to
|
||||
`?`, so transitions to/from the empty state are never scored. But the
|
||||
calibrated 16-state cost matrix (`obikindex::pairwise_cost_matrix`) prices
|
||||
*every* cardinality change — not just `X → ∅` — via a shared
|
||||
`P_cardinality(|A|→|B|)` factor (the same constant this project has called
|
||||
`c_ctx` throughout: "used everywhere a member is gained or lost", see "A
|
||||
concrete Sankoff cost matrix" above). So a transition between two
|
||||
*non-empty* states that differ in how many family members are detected
|
||||
(e.g. `{A} -> {A,C}`, gaining a second sibling) was still carrying the full
|
||||
calibrated gain/loss penalty under `--free-loss`, even though the same
|
||||
low/incomplete-coverage argument for not trusting whole-family absence
|
||||
applies identically one level down: whether a genome shows 1 vs 2 detected
|
||||
members of a family it does carry is exactly as vulnerable to sampling
|
||||
failure as whether the family was detected at all. Fixed by adding a
|
||||
`free_loss: bool` parameter to `pairwise_cost_matrix` that drops the
|
||||
`P_cardinality` factor from the log-probability sum entirely when set,
|
||||
leaving cost driven only by composition matching (shared-base retention,
|
||||
paired substitutions) — never by cardinality difference alone. Unit-tested
|
||||
(`free_loss_ignores_cardinality_transition_probs`): two wildly different
|
||||
cardinality models produce an identical cost matrix once `free_loss` is
|
||||
true. Not yet re-run against a real index to see how much this changes the
|
||||
resulting tree beyond the unit-level guarantee.
|
||||
|
||||
**`+ASC` misbehaves on closely-related taxa, separately from either bug
|
||||
above.** Restricting to the 3-genus, 11-genome subset and running `--iqtree`
|
||||
with `+ASC` converges suspiciously fast compared to the same run without
|
||||
`+ASC` (which takes the expected ~100 NNI iterations, ~1m38s, converges to a
|
||||
stable score, total tree length 7.147, and gives the clean 3-clade
|
||||
topology). Not yet root-caused with a real `+ASC` log in hand (the one run
|
||||
overwrote its own log under the shared `--prefix`, see the operational note
|
||||
below) — working hypothesis, not confirmed: `+ASC`'s correction divides by
|
||||
`P(site variable | current model, current branch lengths)`, recomputed
|
||||
every iteration; on closely-related strains (short branches), a plausible
|
||||
model predicts most sites *should* be constant, so `P(variable)` can be
|
||||
small enough that the correction becomes numerically unstable. Needs a
|
||||
dedicated run (own `--prefix`, log kept) before this is more than a
|
||||
hypothesis.
|
||||
|
||||
**Operational note, not a code bug: IQ-TREE's `--undo`/checkpoint reuse
|
||||
across different alignments crashes.** Reusing `iqtree.ckp.gz` (and the tree
|
||||
it stores) from a run on one taxon set against a regenerated `iqtree.fasta`
|
||||
with a different taxon count crashes with `Assertion 'i >= 0 && i <
|
||||
(int)seq_names.size()' failed` in `Alignment::getSeqName` /
|
||||
`PhyloTree::assignLeafNames` — `--undo` tries to graft a stored tree with
|
||||
the old leaf count onto the new, smaller alignment. Not an `obikmer` bug;
|
||||
avoid by using a fresh `--prefix` (or `-redo`) whenever the underlying
|
||||
alignment's taxon set changes, never `--undo` across them.
|
||||
|
||||
## References
|
||||
|
||||
The Mash mutation-rate model this discussion contrasts with:
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
(Escherichia_coli--CFT073:0.7973338454,(Escherichia_coli--EDL933:0.5132866305,(Escherichia_coli--K-12_MG1655:0.2990687226,Escherichia_coli--K-12_W3110:0.3006288477)100:0.1655629834)69:0.0000007189,(((Klebsiella_pneumoniae--ATCC_13883:0.9661675519,(Klebsiella_pneumoniae--MGH_78578:0.3148953385,Yersinia_ruckeri--YRB:0.3576917943)84:0.4568464394)83:0.0727808943,(Klebsiella_pneumoniae--HS11286:0.1397401144,Proteus_mirabilis--HI4320:0.2192824480)77:0.3686486712)76:0.1162295711,((Salmonella_enterica--AKU_12601:0.6158730439,(Salmonella_enterica--LT2:0.4371443730,Salmonella_enterica--P125109:0.5815549929)100:0.1012287522)100:0.0714169372,Salmonella_enterica--CT18:0.2058143021)76:0.3641489914)100:0.2900565429);
|
||||
@@ -0,0 +1,62 @@
|
||||
#nexus
|
||||
|
||||
BEGIN Taxa;
|
||||
DIMENSIONS ntax=13;
|
||||
TAXLABELS
|
||||
[1] 'Escherichia_coli--CFT073'
|
||||
[2] 'Escherichia_coli--EDL933'
|
||||
[3] 'Escherichia_coli--K-12_MG1655'
|
||||
[4] 'Escherichia_coli--K-12_W3110'
|
||||
[5] 'Klebsiella_pneumoniae--ATCC_13883'
|
||||
[6] 'Klebsiella_pneumoniae--HS11286'
|
||||
[7] 'Klebsiella_pneumoniae--MGH_78578'
|
||||
[8] 'Proteus_mirabilis--HI4320'
|
||||
[9] 'Salmonella_enterica--AKU_12601'
|
||||
[10] 'Salmonella_enterica--CT18'
|
||||
[11] 'Salmonella_enterica--LT2'
|
||||
[12] 'Salmonella_enterica--P125109'
|
||||
[13] 'Yersinia_ruckeri--YRB'
|
||||
;
|
||||
END; [Taxa]
|
||||
|
||||
BEGIN Splits;
|
||||
DIMENSIONS ntax=13 nsplits=34;
|
||||
FORMAT labels=no weights=yes confidences=no intervals=no;
|
||||
MATRIX
|
||||
100 2,
|
||||
100 3,
|
||||
100 4,
|
||||
100 3 4,
|
||||
69 2 3 4,
|
||||
100 5,
|
||||
100 7,
|
||||
16 5 7,
|
||||
100 6,
|
||||
100 13,
|
||||
16 6 13,
|
||||
23 5 6 7 13,
|
||||
100 8,
|
||||
100 10,
|
||||
23 8 10,
|
||||
100 9,
|
||||
100 11,
|
||||
100 12,
|
||||
100 11 12,
|
||||
100 9 11 12,
|
||||
23 8 9 10 11 12,
|
||||
100 1 2 3 4,
|
||||
100 1,
|
||||
84 7 13,
|
||||
83 5 7 13,
|
||||
77 6 8,
|
||||
76 5 6 7 8 13,
|
||||
76 9 10 11 12,
|
||||
31 1 2,
|
||||
1 5 6 7 8,
|
||||
1 10 13,
|
||||
1 9 10 11 12 13,
|
||||
1 5 6 13,
|
||||
1 5 6 8,
|
||||
;
|
||||
END; [Splits]
|
||||
|
||||
+1000
File diff suppressed because it is too large
Load Diff
@@ -46,6 +46,7 @@ nav:
|
||||
- dump: usage/dump.md
|
||||
- annotate: usage/annotate.md
|
||||
- phylo: usage/phylo.md
|
||||
- name-tree: usage/name-tree.md
|
||||
- unitig: usage/unitig.md
|
||||
- estimate: usage/estimate.md
|
||||
- reindex: usage/reindex.md
|
||||
|
||||
+74
-74
@@ -2,166 +2,166 @@ ratio_ceiling: 0.5
|
||||
cardinality_transitions:
|
||||
- from: 0
|
||||
to: 0
|
||||
count: 37622567
|
||||
probability: 0.7765461217831444
|
||||
count: 54999793
|
||||
probability: 0.7934021725308564
|
||||
- from: 0
|
||||
to: 1
|
||||
count: 6977400
|
||||
probability: 0.14401656617768033
|
||||
count: 8873859
|
||||
probability: 0.12801028195383377
|
||||
- from: 0
|
||||
to: 2
|
||||
count: 3744558
|
||||
probability: 0.0772893033240408
|
||||
count: 5305545
|
||||
probability: 0.07653539585976665
|
||||
- from: 0
|
||||
to: 3
|
||||
count: 91083
|
||||
probability: 0.0018799926759482985
|
||||
count: 126109
|
||||
probability: 0.001819191475424167
|
||||
- from: 0
|
||||
to: 4
|
||||
count: 12985
|
||||
probability: 0.0002680160391861122
|
||||
count: 16149
|
||||
probability: 0.00023295818011898336
|
||||
- from: 1
|
||||
to: 0
|
||||
count: 6977400
|
||||
probability: 0.820630384894398
|
||||
count: 8873859
|
||||
probability: 0.8510749571674425
|
||||
- from: 1
|
||||
to: 1
|
||||
count: 1018042
|
||||
probability: 0.1197346000370715
|
||||
count: 1041044
|
||||
probability: 0.09984455215137214
|
||||
- from: 1
|
||||
to: 2
|
||||
count: 502624
|
||||
probability: 0.05911493200578466
|
||||
count: 507272
|
||||
probability: 0.048651493749477304
|
||||
- from: 1
|
||||
to: 3
|
||||
count: 4197
|
||||
probability: 0.000493620220340211
|
||||
count: 4248
|
||||
probability: 0.00040741760918753563
|
||||
- from: 1
|
||||
to: 4
|
||||
count: 225
|
||||
probability: 0.000026462842405658202
|
||||
probability: 0.00002157932252052625
|
||||
- from: 2
|
||||
to: 0
|
||||
count: 3744558
|
||||
probability: 0.8804244781723005
|
||||
count: 5305545
|
||||
probability: 0.9117891245815953
|
||||
- from: 2
|
||||
to: 1
|
||||
count: 502624
|
||||
probability: 0.11817749195415704
|
||||
count: 507272
|
||||
probability: 0.08717767784549091
|
||||
- from: 2
|
||||
to: 2
|
||||
count: 5162
|
||||
probability: 0.0012136949558066439
|
||||
count: 5227
|
||||
probability: 0.0008982907041949506
|
||||
- from: 2
|
||||
to: 3
|
||||
count: 688
|
||||
probability: 0.0001617632951559417
|
||||
count: 689
|
||||
probability: 0.00011840870388182915
|
||||
- from: 2
|
||||
to: 4
|
||||
count: 96
|
||||
probability: 0.00002257162257989884
|
||||
probability: 0.000016498164836945715
|
||||
- from: 3
|
||||
to: 0
|
||||
count: 91083
|
||||
probability: 0.9478728718311618
|
||||
count: 126109
|
||||
probability: 0.9614164824273843
|
||||
- from: 3
|
||||
to: 1
|
||||
count: 4197
|
||||
probability: 0.04367689297756317
|
||||
count: 4248
|
||||
probability: 0.03238545399100404
|
||||
- from: 3
|
||||
to: 2
|
||||
count: 688
|
||||
probability: 0.007159805186696083
|
||||
count: 689
|
||||
probability: 0.005252725470763132
|
||||
- from: 3
|
||||
to: 3
|
||||
count: 98
|
||||
probability: 0.0010198559713607792
|
||||
probability: 0.0007471220553480217
|
||||
- from: 3
|
||||
to: 4
|
||||
count: 26
|
||||
probability: 0.00027057403321816594
|
||||
probability: 0.00019821605550049553
|
||||
- from: 4
|
||||
to: 0
|
||||
count: 12985
|
||||
probability: 0.972950696838004
|
||||
count: 16149
|
||||
probability: 0.9781344639612356
|
||||
- from: 4
|
||||
to: 1
|
||||
count: 225
|
||||
probability: 0.01685898396523303
|
||||
probability: 0.013628104179285281
|
||||
- from: 4
|
||||
to: 2
|
||||
count: 96
|
||||
probability: 0.0071931664918327585
|
||||
probability: 0.00581465778316172
|
||||
- from: 4
|
||||
to: 3
|
||||
count: 26
|
||||
probability: 0.0019481492582047056
|
||||
probability: 0.0015748031496062992
|
||||
- from: 4
|
||||
to: 4
|
||||
count: 14
|
||||
probability: 0.0010490034467256107
|
||||
probability: 0.0008479709267110841
|
||||
composition_transitions:
|
||||
- from: 'A'
|
||||
to: 'A'
|
||||
count: 157182
|
||||
probability: 0.46648741052031767
|
||||
count: 161774
|
||||
probability: 0.4701762136303263
|
||||
- from: 'A'
|
||||
to: 'C'
|
||||
count: 27305
|
||||
probability: 0.08103624298111281
|
||||
count: 27788
|
||||
probability: 0.0807624007835592
|
||||
- from: 'A'
|
||||
to: 'G'
|
||||
count: 128548
|
||||
probability: 0.38150693875612857
|
||||
count: 130019
|
||||
probability: 0.3778842157577942
|
||||
- from: 'A'
|
||||
to: 'T'
|
||||
count: 23913
|
||||
probability: 0.07096940774244097
|
||||
count: 24490
|
||||
probability: 0.07117716982832031
|
||||
- from: 'C'
|
||||
to: 'A'
|
||||
count: 27305
|
||||
probability: 0.07823849488677553
|
||||
count: 27788
|
||||
probability: 0.0780165140757087
|
||||
- from: 'C'
|
||||
to: 'C'
|
||||
count: 179118
|
||||
probability: 0.5132365034656459
|
||||
count: 184047
|
||||
probability: 0.5167232390273484
|
||||
- from: 'C'
|
||||
to: 'G'
|
||||
count: 19392
|
||||
probability: 0.055564947549692406
|
||||
count: 19622
|
||||
probability: 0.05508996830263265
|
||||
- from: 'C'
|
||||
to: 'T'
|
||||
count: 123182
|
||||
probability: 0.3529600540978862
|
||||
count: 124724
|
||||
probability: 0.3501702785943102
|
||||
- from: 'G'
|
||||
to: 'A'
|
||||
count: 128548
|
||||
probability: 0.3637269077214298
|
||||
count: 130019
|
||||
probability: 0.3610094570655886
|
||||
- from: 'G'
|
||||
to: 'C'
|
||||
count: 19392
|
||||
probability: 0.05486971555009776
|
||||
count: 19622
|
||||
probability: 0.05448224926003876
|
||||
- from: 'G'
|
||||
to: 'G'
|
||||
count: 179245
|
||||
probability: 0.5071742039901647
|
||||
count: 183951
|
||||
probability: 0.5107565097152884
|
||||
- from: 'G'
|
||||
to: 'T'
|
||||
count: 26234
|
||||
probability: 0.07422917273830779
|
||||
count: 26562
|
||||
probability: 0.07375178395908417
|
||||
- from: 'T'
|
||||
to: 'A'
|
||||
count: 23913
|
||||
probability: 0.07307212790143376
|
||||
count: 24490
|
||||
probability: 0.07335783586895636
|
||||
- from: 'T'
|
||||
to: 'C'
|
||||
count: 123182
|
||||
probability: 0.3764132839524281
|
||||
count: 124724
|
||||
probability: 0.37360076443118473
|
||||
- from: 'T'
|
||||
to: 'G'
|
||||
count: 26234
|
||||
probability: 0.08016452153080807
|
||||
count: 26562
|
||||
probability: 0.07956434611479049
|
||||
- from: 'T'
|
||||
to: 'T'
|
||||
count: 153923
|
||||
probability: 0.4703500666153301
|
||||
count: 158067
|
||||
probability: 0.4734770535850684
|
||||
|
||||
@@ -129,7 +129,23 @@ fn best_pairing_cost(lost: &[u8], gained: &[u8], p_comp: &[[f64; 4]; 4]) -> f64
|
||||
/// symmetrised (`(cost(A,B)+cost(B,A))/2` — the row-normalised `P` is not
|
||||
/// symmetric in general, but a Sankoff parsimony cost must be, so the
|
||||
/// score is independent of where an unrooted tree gets rooted).
|
||||
pub fn pairwise_cost_matrix(p_card: &[[f64; 5]; 5], p_comp: &[[f64; 4]; 4]) -> [[f64; 16]; 16] {
|
||||
///
|
||||
/// `free_loss`: drop the `P_cardinality(|A|→|B|)` factor entirely (never
|
||||
/// added to `log_p`) — the same low/incomplete-coverage argument that
|
||||
/// justifies recoding whole-family non-detection as `?` (see
|
||||
/// `docmd/theory/evolutionary_distances.md`, "Locus dropout under
|
||||
/// incomplete coverage") applies one level down: whether a genome shows 1
|
||||
/// vs 2 (etc.) detected members of a *present* family is exactly as
|
||||
/// vulnerable to sampling failure as whether the family was detected at
|
||||
/// all. Without this, `∅`-involving transitions are neutralised (via the
|
||||
/// `?` recoding, bypassing this matrix's row/column 0 entirely) but
|
||||
/// cardinality changes *between two otherwise-detected, non-empty* states
|
||||
/// (e.g. `{A} -> {A,C}`) still carried the same calibrated
|
||||
/// `P_cardinality` penalty as any other gain/loss — inconsistent with
|
||||
/// `--free-loss`'s own rationale. With the factor dropped, cost is driven
|
||||
/// only by composition matching (shared-base retention and paired
|
||||
/// substitutions), never by a state pair's cardinality difference alone.
|
||||
pub fn pairwise_cost_matrix(p_card: &[[f64; 5]; 5], p_comp: &[[f64; 4]; 4], free_loss: bool) -> [[f64; 16]; 16] {
|
||||
let mut raw = [[0.0f64; 16]; 16];
|
||||
for a in 0u8..16 {
|
||||
for b in 0u8..16 {
|
||||
@@ -138,10 +154,12 @@ pub fn pairwise_cost_matrix(p_card: &[[f64; 5]; 5], p_comp: &[[f64; 4]; 4]) -> [
|
||||
let gained: Vec<u8> = (0..4).filter(|&i| b & (1 << i) != 0 && a & (1 << i) == 0).collect();
|
||||
|
||||
let mut log_p = 0.0; // accumulate ln(P), so 0.0 = probability 1
|
||||
let card_a = a.count_ones() as usize;
|
||||
let card_b = b.count_ones() as usize;
|
||||
let p_c = p_card[card_a][card_b];
|
||||
log_p += if p_c > 0.0 { p_c.ln() } else { f64::NEG_INFINITY };
|
||||
if !free_loss {
|
||||
let card_a = a.count_ones() as usize;
|
||||
let card_b = b.count_ones() as usize;
|
||||
let p_c = p_card[card_a][card_b];
|
||||
log_p += if p_c > 0.0 { p_c.ln() } else { f64::NEG_INFINITY };
|
||||
}
|
||||
|
||||
for i in 0..4u8 {
|
||||
if shared & (1 << i) != 0 {
|
||||
@@ -231,7 +249,7 @@ mod tests {
|
||||
[0.15, 0.05, 0.7, 0.1],
|
||||
[0.05, 0.15, 0.1, 0.7],
|
||||
];
|
||||
let cost = pairwise_cost_matrix(&p_card, &p_comp);
|
||||
let cost = pairwise_cost_matrix(&p_card, &p_comp, false);
|
||||
for a in 0..16 {
|
||||
assert_eq!(cost[a][a], 0.0);
|
||||
for b in 0..16 {
|
||||
@@ -239,4 +257,36 @@ mod tests {
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn free_loss_ignores_cardinality_transition_probs() {
|
||||
// A skewed cardinality model (cardinality change made artificially
|
||||
// expensive) must have zero effect on the cost matrix once
|
||||
// `free_loss` is set — the whole point of the flag.
|
||||
let p_card_uniform = [[0.2f64; 5]; 5];
|
||||
let p_card_skewed = [
|
||||
[0.96, 0.01, 0.01, 0.01, 0.01],
|
||||
[0.01, 0.96, 0.01, 0.01, 0.01],
|
||||
[0.01, 0.01, 0.96, 0.01, 0.01],
|
||||
[0.01, 0.01, 0.01, 0.96, 0.01],
|
||||
[0.01, 0.01, 0.01, 0.01, 0.96],
|
||||
];
|
||||
let p_comp = [
|
||||
[0.7, 0.1, 0.15, 0.05],
|
||||
[0.1, 0.7, 0.05, 0.15],
|
||||
[0.15, 0.05, 0.7, 0.1],
|
||||
[0.05, 0.15, 0.1, 0.7],
|
||||
];
|
||||
let cost_uniform = pairwise_cost_matrix(&p_card_uniform, &p_comp, true);
|
||||
let cost_skewed = pairwise_cost_matrix(&p_card_skewed, &p_comp, true);
|
||||
for a in 0..16 {
|
||||
for b in 0..16 {
|
||||
assert!(
|
||||
(cost_uniform[a][b] - cost_skewed[a][b]).abs() < 1e-9,
|
||||
"cost[{a}][{b}] differs between cardinality models under free_loss: {} vs {}",
|
||||
cost_uniform[a][b], cost_skewed[a][b],
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -9,6 +9,7 @@ pub mod dump;
|
||||
pub mod estimate;
|
||||
pub mod index;
|
||||
pub mod merge;
|
||||
pub mod nametree;
|
||||
pub mod query;
|
||||
pub mod reindex;
|
||||
pub mod superkmer;
|
||||
|
||||
@@ -0,0 +1,146 @@
|
||||
use std::path::{Path, PathBuf};
|
||||
|
||||
use clap::Args;
|
||||
use tracing::info;
|
||||
|
||||
// ── Translate a numerically-labelled tree export back to real taxon names ──
|
||||
//
|
||||
// TNT/PhyG write bare numeric leaf labels (1-based, in the same order as the
|
||||
// FASTA fed to them) — this reads that order back from the FASTA header
|
||||
// line and emits a NEXUS `translate` table alongside the tree(s), unchanged
|
||||
// otherwise. Readable directly by FigTree/PearTree/`ape` etc.
|
||||
|
||||
#[derive(Args)]
|
||||
pub struct NameTreeArgs {
|
||||
/// Tree file to translate: a TNT-style NEXUS export (`tree NAME = [&U]
|
||||
/// ...;`, topology on the same or the next line) or a plain Newick file
|
||||
/// (single `(...);` tree, no header)
|
||||
pub tree: PathBuf,
|
||||
|
||||
/// FASTA file whose record order gives the numeric taxon labels
|
||||
/// (1-based) — typically the `_sankoff.fasta`/`_snp.fasta` used to
|
||||
/// produce `tree`
|
||||
#[arg(long)]
|
||||
pub fasta: PathBuf,
|
||||
|
||||
/// Output NEXUS file (taxa block + translate table + tree(s), topology
|
||||
/// unchanged)
|
||||
#[arg(short, long)]
|
||||
pub output: PathBuf,
|
||||
}
|
||||
|
||||
pub fn run(args: NameTreeArgs) {
|
||||
let labels = read_fasta_labels(&args.fasta);
|
||||
if labels.is_empty() {
|
||||
eprintln!("error: no FASTA headers found in {}", args.fasta.display());
|
||||
std::process::exit(1);
|
||||
}
|
||||
|
||||
let content = std::fs::read_to_string(&args.tree).unwrap_or_else(|e| {
|
||||
eprintln!("error reading {}: {e}", args.tree.display());
|
||||
std::process::exit(1);
|
||||
});
|
||||
let trees = extract_trees(&content);
|
||||
if trees.is_empty() {
|
||||
eprintln!("error: no tree found in {}", args.tree.display());
|
||||
std::process::exit(1);
|
||||
}
|
||||
|
||||
write_named_nexus(&labels, &trees, &args.output);
|
||||
}
|
||||
|
||||
fn read_fasta_labels(path: &Path) -> Vec<String> {
|
||||
let content = std::fs::read_to_string(path).unwrap_or_else(|e| {
|
||||
eprintln!("error reading {}: {e}", path.display());
|
||||
std::process::exit(1);
|
||||
});
|
||||
content.lines()
|
||||
.filter(|l| l.starts_with('>'))
|
||||
.map(|l| {
|
||||
let header = &l[1..];
|
||||
// `obifastwrite::write_record` appends a ` {json}` annotation —
|
||||
// not part of the taxon name.
|
||||
match header.find(" {") {
|
||||
Some(pos) => header[..pos].to_string(),
|
||||
None => header.to_string(),
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Finds every `tree NAME = [&U] TOPOLOGY;` (rooting comment optional,
|
||||
/// topology on the same line or the next non-empty one), or — if none of
|
||||
/// that syntax is found — treats the whole file as one bare Newick tree.
|
||||
fn extract_trees(content: &str) -> Vec<(String, String)> {
|
||||
let lines: Vec<&str> = content.lines().collect();
|
||||
let mut trees = Vec::new();
|
||||
let mut i = 0;
|
||||
while i < lines.len() {
|
||||
let line = lines[i].trim();
|
||||
if let Some(rest) = line.strip_prefix("tree ") {
|
||||
if let Some(eq_pos) = rest.find('=') {
|
||||
let name = rest[..eq_pos].trim().to_string();
|
||||
let mut after_eq = rest[eq_pos + 1..].trim();
|
||||
if after_eq.starts_with('[') {
|
||||
if let Some(close) = after_eq.find(']') {
|
||||
after_eq = after_eq[close + 1..].trim();
|
||||
}
|
||||
}
|
||||
let topo = if after_eq.starts_with('(') {
|
||||
after_eq.to_string()
|
||||
} else {
|
||||
i += 1;
|
||||
while i < lines.len() && lines[i].trim().is_empty() {
|
||||
i += 1;
|
||||
}
|
||||
lines.get(i).map(|s| s.trim().to_string()).unwrap_or_default()
|
||||
};
|
||||
if topo.starts_with('(') {
|
||||
trees.push((name, topo));
|
||||
}
|
||||
}
|
||||
}
|
||||
i += 1;
|
||||
}
|
||||
|
||||
if trees.is_empty() {
|
||||
let trimmed = content.trim();
|
||||
if trimmed.starts_with('(') && trimmed.ends_with(';') {
|
||||
trees.push(("tree_1".to_string(), trimmed.to_string()));
|
||||
}
|
||||
}
|
||||
|
||||
trees
|
||||
}
|
||||
|
||||
fn write_named_nexus(labels: &[String], trees: &[(String, String)], output: &Path) {
|
||||
let mut out = String::new();
|
||||
out.push_str("#NEXUS\n\n");
|
||||
out.push_str("begin taxa;\n");
|
||||
out.push_str(&format!(" dimensions ntax={};\n", labels.len()));
|
||||
out.push_str(" taxlabels\n");
|
||||
for lab in labels {
|
||||
out.push_str(&format!(" {lab}\n"));
|
||||
}
|
||||
out.push_str(" ;\nend;\n\n");
|
||||
out.push_str("begin trees;\n");
|
||||
out.push_str(" translate\n");
|
||||
let tr_lines: Vec<String> = labels.iter().enumerate()
|
||||
.map(|(i, lab)| format!(" {} {lab}", i + 1))
|
||||
.collect();
|
||||
out.push_str(&tr_lines.join(",\n"));
|
||||
out.push_str(";\n");
|
||||
for (name, topo) in trees {
|
||||
out.push_str(&format!(" tree {name} = [&U] {topo}\n"));
|
||||
}
|
||||
out.push_str("end;\n");
|
||||
|
||||
std::fs::write(output, &out).unwrap_or_else(|e| {
|
||||
eprintln!("error writing {}: {e}", output.display());
|
||||
std::process::exit(1);
|
||||
});
|
||||
info!(
|
||||
"named tree(s) → {} ({} tree{}, {} taxa)",
|
||||
output.display(), trees.len(), if trees.len() == 1 { "" } else { "s" }, labels.len(),
|
||||
);
|
||||
}
|
||||
@@ -184,7 +184,7 @@ pub fn run(args: PhyloArgs) {
|
||||
});
|
||||
let p_card = cardinality_transition_probs(&card_tally);
|
||||
let p_comp = composition_transition_probs(&base_tally);
|
||||
let matrix = pairwise_cost_matrix(&p_card, &p_comp);
|
||||
let matrix = pairwise_cost_matrix(&p_card, &p_comp, args.free_loss);
|
||||
write_sankoff_matrix_csv(&matrix, &args.output);
|
||||
write_sankoff_params(&card_tally, &p_card, &base_tally, &p_comp, args.sankoff_ratio_ceiling, &args.output);
|
||||
|
||||
|
||||
@@ -33,6 +33,9 @@ enum Commands {
|
||||
/// Compute pairwise evolutionary-distance proxies between genomes (metric matrix, NJ/UPGMA
|
||||
/// trees, SNP/Sankoff calibration, TNT/PhyG/IQ-TREE exports)
|
||||
Phylo(cmd::phylo::PhyloArgs),
|
||||
/// Translate a numerically-labelled tree export (TNT/PhyG) back to real taxon names, from
|
||||
/// the FASTA that produced it
|
||||
NameTree(cmd::nametree::NameTreeArgs),
|
||||
/// Dump unitigs from a built index to stdout (debug)
|
||||
Unitig(cmd::unitig::UnitigArgs),
|
||||
/// Estimate approximate-index parameters (z, evidence bits, FP rates) before indexing
|
||||
@@ -73,6 +76,7 @@ fn main() {
|
||||
Commands::Query(args) => cmd::query::run(args),
|
||||
Commands::Annotate(args) => cmd::annotate::run(args),
|
||||
Commands::Phylo(args) => cmd::phylo::run(args),
|
||||
Commands::NameTree(args) => cmd::nametree::run(args),
|
||||
Commands::Unitig(args) => cmd::unitig::run(args),
|
||||
Commands::Estimate(args) => cmd::estimate::run(args),
|
||||
Commands::Reindex(args) => cmd::reindex::run(args),
|
||||
|
||||
Reference in New Issue
Block a user