obikmer/UserDocMD/usage/phylo.md

# phylo

Compute pairwise evolutionary-distance proxies between the genomes stored in an index — a plain distance matrix, optionally trees (NJ/UPGMA), and optionally a central-position SNP model with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).

```bash
obikmer phylo INDEX [OPTIONS]
```

## Arguments

| Argument | Description |
|---|---|
| `INDEX` | Index directory |

## Distance matrix

| Option | Default | Description |
|---|---|---|
| `--metric` | `jaccard` | One of `jaccard`, `mash`, `hamming`, `bray-curtis`, `relfreq-bray-curtis`, `euclidean`, `relfreq-euclidean`, `hellinger`, `hellinger-euclidean` |
| `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for Jaccard/Mash on a count index |
| `--shared-kmers` | off | Also write the shared-kmer count matrix |
| `--nj` | off | Compute and write a Neighbor-Joining tree (Newick) |
| `--upgma` | off | Compute and write a UPGMA tree (Newick) |
| `-o, --output` | none (stdout) | Output file prefix; without it, the distance matrix is printed to stdout as CSV |

`hamming` requires a presence/absence index. All other metrics work on either index type; on a presence index, `jaccard`/`mash`/`hamming` are the only ones available.

### Metric definitions

- **jaccard**: $D = 1 - \dfrac{|A \cap B|}{|A \cup B|}$ over the sets of kmers present in each genome.
- **mash**: derived from the Jaccard distance via $D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)$ where $J = 1 - D_{\text{jaccard}}$ and $k$ is the index's kmer size; clamped to 1.0 when $J \le 0$.
- **hamming**: number of kmer positions where presence differs between the two genomes (presence index only, not normalized): $D = \sum_i \mathbb{1}[a_i \ne b_i]$.
- **bray-curtis**: $D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$ on raw per-kmer counts.
- **relfreq-bray-curtis**: the same formula computed on per-genome relative frequencies $p_i = c_i / \sum_j c_j$ instead of raw counts.
- **euclidean**: $D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$ on raw counts.
- **relfreq-euclidean**: the same formula on relative frequencies.
- **hellinger**: $D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$ on relative frequencies, bounded in $[0, 1]$.
- **hellinger-euclidean**: the unnormalized variant, $D = \sqrt{2} \times D_{\text{hellinger}}$.

## Central-position SNP model

This is a separate operation from the distance-matrix computation above: if any option below is used, no `--metric` matrix is computed in the same invocation.

A **family** is the set of up to 4 kmers that share identical flanking sequence and differ only at the exact central base. Because $k$ is odd, the central position is well defined and maps to itself under reverse complementation. All computations below first require building the **sibling annex**, an index-wide record of which of the 4 possible central bases are observed at each family, across every genome.

| Option | Default | Description |
|---|---|---|
| `--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

A family is eligible for a genome pair $(i, j)$ only if genome $i$ carries exactly one of the family's observed forms (single-copy, unambiguous) and genome $j$ also carries exactly one. A genome carrying more than one form at a locus makes that locus ineligible for any pair involving it.

`--raw-snp-distance` tallies, over every eligible locus of every genome pair, $n_{\text{snp}}$ (the two genomes' single forms differ) versus $n_{\text{shared}}$ (they agree — this includes invariant families). The output ratio is $\hat{p} = \dfrac{n_{\text{snp}}}{n_{\text{snp}} + n_{\text{shared}}}$.

`--snp` restricts itself to *variable* families (family size $\ge 2$) and writes one FASTA record per genome, one column per family, IUPAC-coded from each genome's presence mask at that family (a single form → the plain base; several forms → the matching IUPAC ambiguity code; no form → `-`).

`--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`) |
| `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands) |

### The 16-state model

Each family is treated as a character with 16 possible states: one per subset of the 4 possible central bases actually observed (including the empty subset). Calibration combines two tallies, both restricted to genome pairs at or below `--sankoff-ratio-ceiling`:

- a $5 \times 5$ transition matrix over family cardinality (0–4 observed forms) between paired genomes, and
- a $4 \times 4$ base-substitution transition matrix from unambiguous single-copy loci,

which are combined into a row-normalized $16 \times 16$ transition probability matrix $P$, converted to a symmetric cost matrix via $\text{cost}(a,b) = -\ln P(a,b)$.

`--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. `--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

All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment, recoded for the target tool:

- **`--tnt`**: a self-contained TNT script (alignment recoded to TNT's fixed 16-symbol alphabet, integer-scaled cost matrix re-closed to a metric, a default search block).
- **`--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 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
```