From 3da501349b49350995bae1c8a658ca15c3db74dd Mon Sep 17 00:00:00 2001
From: Eric Coissac <eric@coissac.eu>
Date: Fri, 14 Aug 2026 07:43:07 +0200
Subject: [PATCH] docs: clarify phylogenetic output files and tool usage

Restructure the output files section into categorized subsections with tables. Add explicit mappings between command-line options and generated files. Define CSV matrix conventions, clarify mathematical formulas for distance calculations, and document execution commands for external phylogenetic tools.
---
 UserDocMD/usage/phylo.md | 84 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 83 insertions(+), 1 deletion(-)

diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md
index 378a33d1..2d6b5e88 100644
--- a/UserDocMD/usage/phylo.md
+++ b/UserDocMD/usage/phylo.md
@@ -107,4 +107,86 @@ TNT and PhyG both write trees with bare numeric leaf labels (`1`, `2`, …, in t
 
 ## 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.
+With `-o/--output PREFIX`, the relevant subset of the files below is written. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout. All matrices use genome labels (from the index metadata) as row/column headers, in index order; all CSVs are comma-separated with a header row.
+
+### Distance matrix
+
+| File | Written by | Format | Content |
+|---|---|---|---|
+| `<prefix>_dist.csv` | always | CSV matrix | the `--metric` distance, 6 decimals, symmetric, diagonal 0 |
+| `<prefix>_shared.csv` | `--shared-kmers` | CSV matrix | shared-kmer count per genome pair (integers) |
+| `<prefix>_nj.nwk` | `--nj` | Newick | Neighbor-Joining tree, branch lengths from the `--metric` matrix |
+| `<prefix>_upgma.nwk` | `--upgma` | Newick | UPGMA tree, same matrix |
+
+Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below): header `genome,<label1>,<label2>,...`, one data row per genome, `<label>,<value1>,<value2>,...`.
+
+### Central-position SNP model
+
+| File | Written by | Format | Content |
+|---|---|---|---|
+| `<prefix>_siblings.csv` | `--sibling-stats` | CSV table | family-size distribution, per genome and global |
+| `<prefix>_rawsnp.csv` | `--raw-snp-distance` | CSV matrix | single-copy central-SNP p-distance ($\hat p$), or `NA` |
+| `<prefix>_rawsnp_counts.csv` | `--raw-snp-counts` | CSV table | per-pair diagnostic counts behind `_rawsnp.csv` |
+| `<prefix>_snp.fasta` | `--snp` | FASTA | SNP-only pseudo-alignment, IUPAC-coded |
+| `<prefix>_family_overlap.csv` | `--family-overlap` | CSV matrix | variable families both genomes of a pair carry a call for |
+
+**`_siblings.csv`** — family size = number of distinct central bases observed at a family (1–4), not "sibling count" (0–3).
+
+| Column | Meaning |
+|---|---|
+| `genome` | genome label, or the literal `global` for the last row |
+| `1`, `2`, `3`, `4` | for a genome row: number of families of that size where the genome carries ≥ 1 member. For the `global` row: the actual deduplicated family-size histogram — **not** the sum of the rows above (a family shared by several genomes would otherwise be counted once per genome) |
+
+**`_rawsnp.csv`** — same matrix layout as `_dist.csv`; each cell is $\hat p = n_{\text{snp}}/(n_{\text{snp}}+n_{\text{shared}})$, 6 decimals, or `NA` when the pair has zero eligible loci (distinguishes "identical everywhere eligible" from "nothing eligible at all").
+
+**`_rawsnp_counts.csv`** — one row per unordered genome pair (not a matrix), the counts `_rawsnp.csv`'s ratio is computed from:
+
+| Column | Meaning |
+|---|---|
+| `genome_a`, `genome_b` | the pair |
+| `n_snp` | eligible loci where the two genomes' single forms differ |
+| `n_shared` | eligible loci where they agree (includes invariant families) |
+| `n_eligible` | `n_snp + n_shared` |
+| `ratio` | $\hat p$ = `n_snp / n_eligible`, or `NA` if `n_eligible = 0` |
+
+**`_snp.fasta`** — one record per non-excluded genome, one column per variable family (family size ≥ 2), header carries an `n_sites` annotation. Each site is IUPAC-coded from the genome's presence mask at that family: single observed form → plain base; several forms → matching IUPAC ambiguity code; no form → `-`.
+
+**`_family_overlap.csv`** — same matrix layout as `_dist.csv`; cell `[i][j]` = number of `_snp.fasta` columns where both genome `i` and `j` carry a call (neither is `-`). Diagonal `[i][i]` is kept (not skipped): it holds genome `i`'s own total variable-family count.
+
+### Sankoff calibration and exports
+
+| File | Written by | Format | Content |
+|---|---|---|---|
+| `<prefix>_sankoff_matrix.csv` | `--sankoff`/`--tnt`/`--phyg`/`--iqtree` | CSV matrix | calibrated 16×16 cost matrix |
+| `<prefix>_sankoff_params.yaml` | same flags | YAML | calibration report (raw tallies + derived probabilities) |
+| `<prefix>_sankoff.fasta` | same flags | FASTA | Sankoff-recoded pseudo-alignment |
+| `<prefix>_sankoff.tnt` | `--tnt` | TNT script | ready-to-run parsimony search |
+| `<prefix>_sankoff.tcm` | `--phyg` | PhyG TCM | cost matrix in PhyG's own format |
+| `<prefix>_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search |
+| `<prefix>_iqtree.model` | `--iqtree` | IQ-TREE model file | custom ML substitution model |
+| `<prefix>_iqtree.fasta` | `--iqtree` | FASTA | alignment recoded for that model |
+
+**`_sankoff_matrix.csv`** — header `state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N`: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), `0` standing for the empty/absent state (not `-`, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost $-\ln P(a,b)$, 4 decimals.
+
+**`_sankoff_params.yaml`** — everything the calibration estimated, structured so it can be reloaded rather than re-parsed:
+
+| Key | Meaning |
+|---|---|
+| `ratio_ceiling` | the `--sankoff-ratio-ceiling` value used |
+| `cardinality_transitions` | 5×5 list of `{from, to, count, probability}`, family cardinality (0–4 observed forms) |
+| `composition_transitions` | 4×4 list of `{from, to, count, probability}`, base letters `A/C/G/T`, single-copy substitutions |
+
+**`_sankoff.fasta`** — same sites as `_snp.fasta`, recoded to match `_sankoff_matrix.csv`'s alphabet: absent state is `0` (or `?` under `--free-loss`). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.
+
+**`_sankoff.tnt`** (`--tnt`) — self-contained TNT script: `xread` block (alignment recoded to TNT's fixed `0-9A-F` alphabet), an integer-scaled (`--sankoff-cost-scale`) and metric-closed `smatrix`, a default `hold 20; mult; export` search. Run with `printf 'proc <path>;\nquit;\n' | tnt`. Produces `<prefix>_sankoff.tre` (bare numeric leaf labels, order matching `_sankoff.fasta`) — feed both into [`name-tree`](name-tree.md) to recover taxon names.
+
+**`_sankoff.tcm`** (`--phyg`) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state `0`, since it's never actually triggered).
+
+**`_sankoff.pg`** (`--phyg`) — script: `read(prefasta:..., tcm:...)` against `_sankoff.fasta`/`_sankoff.tcm`, a default 300s/4-instance `search`, `report(...)` writing `<prefix>_sankoff.tre` (bare numeric labels, as for `--tnt`). Run with `phyg` from the output directory (the script uses relative file names). Feed the tree plus `_sankoff.fasta` into [`name-tree`](name-tree.md) for taxon names.
+
+**`_iqtree.model`** (`--iqtree`) — lower-triangular exchangeability matrix $R(a,b) = e^{-\text{cost}(a,b)}$ (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered `0..k-1`.
+
+**`_iqtree.fasta`** (`--iqtree`) — alignment recoded to that same compact `0..k-1` alphabet (symbols `0-9A-F`). Under `--free-loss`, non-detection becomes `?` and columns left non-informative once missing calls are ignored are dropped first (required for `+ASC`). Run with:
+```
+iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
+```