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.
This commit is contained in:
Eric Coissac
@@ -27,3 +27,12 @@ benchmark/specific_index_presence
|
||||
TNT
|
||||
phyg
|
||||
*.tnt
|
||||
*.tre
|
||||
*.phy
|
||||
*.treefile
|
||||
*.bionj
|
||||
*.iqtree
|
||||
*.mldist
|
||||
*.parstree
|
||||
*.ckp.gz
|
||||
*.model
|
||||
|
||||
@@ -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: `<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>_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.
|
||||
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.
|
||||
|
||||
@@ -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:
|
||||
|
||||
@@ -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
|
||||
+74
-74
@@ -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
|
||||
|
||||
@@ -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<MetricArg> 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<String>,
|
||||
|
||||
/// 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<f64>,
|
||||
|
||||
/// 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 (`<prefix>_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 <prefix>_sankoff.tnt, a ready-to-run TNT script (`proc
|
||||
/// <file>;`) 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 <prefix>_sankoff.tcm and <prefix>_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 <prefix>_iqtree.model and <prefix>_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: <prefix>_dist.csv, <prefix>_shared.csv,
|
||||
/// <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,
|
||||
/// <prefix>_nj.nwk, <prefix>_upgma.nwk.
|
||||
/// If omitted, the distance matrix is written to stdout.
|
||||
#[arg(short, long)]
|
||||
pub output: Option<PathBuf>,
|
||||
}
|
||||
@@ -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<Vec<u64>> {
|
||||
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<u64>], 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<PathBuf>) {
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -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<bool> = (0..n_sites).map(|site| {
|
||||
let mut first: Option<u8> = 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<PathBuf>,
|
||||
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<u8> = 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<PathBuf>,
|
||||
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())
|
||||
|
||||
@@ -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<MetricArg> 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<String>,
|
||||
|
||||
/// 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 <prefix>_sankoff.tnt, a ready-to-run TNT script (`proc
|
||||
/// <file>;`) 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 <prefix>_sankoff.tcm and <prefix>_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 <prefix>_iqtree.model and <prefix>_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: <prefix>_dist.csv, <prefix>_shared.csv,
|
||||
/// <prefix>_siblings.csv, <prefix>_rawsnp.csv, <prefix>_rawsnp_counts.csv,
|
||||
/// <prefix>_snp.fasta,
|
||||
/// <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,
|
||||
/// <prefix>_nj.nwk, <prefix>_upgma.nwk.
|
||||
/// If omitted, the distance matrix is written to stdout.
|
||||
#[arg(short, long)]
|
||||
pub output: Option<PathBuf>,
|
||||
}
|
||||
|
||||
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<PathBuf>) {
|
||||
// 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<PathBuf>) {
|
||||
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<PathBuf>) {
|
||||
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<PathBuf>) {
|
||||
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<f64>, 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<String> = names.to_vec();
|
||||
// height of each node: leaves = 0, internal = dissimilarity/2
|
||||
let mut heights: Vec<f64> = 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())
|
||||
}
|
||||
|
||||
@@ -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<PathBuf>) {
|
||||
// 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<PathBuf>) {
|
||||
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<PathBuf>) {
|
||||
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<PathBuf>) {
|
||||
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<f64>, 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<String> = names.to_vec();
|
||||
// height of each node: leaves = 0, internal = dissimilarity/2
|
||||
let mut heights: Vec<f64> = 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())
|
||||
}
|
||||
@@ -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<PathBuf>) {
|
||||
pub(super) fn write_sankoff_alignment_fasta(alignment: &SnpAlignment, labels: &[String], output: &Option<PathBuf>, 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<u8> = seq.iter().map(|&b| if b == b'-' { b'0' } else { b }).collect();
|
||||
let recoded: Vec<u8> = 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);
|
||||
|
||||
@@ -28,6 +28,7 @@ pub(super) fn write_sankoff_tnt(
|
||||
labels: &[String],
|
||||
output: &Option<PathBuf>,
|
||||
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();
|
||||
|
||||
Reference in New Issue
Block a user