# phylo Compute pairwise distances between the genomes stored in an index, optionally build trees (NJ/UPGMA) from them, and optionally calibrate a 16-state parsimony model for a central-position SNP character with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE). ```bash obikmer phylo INDEX [OPTIONS] ``` ## Arguments | Argument | Description | |---|---| | `INDEX` | Index directory | ## Distance matrix (`--distance`) | Option | Default | Description | |---|---|---| | `--distance` | `jaccard` | See the two tables below for the full list of accepted values | | `--gamma-shape ALPHA` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). No effect on the other values; rejected if given together with a value that doesn't support it | | `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for `jaccard`/`mash` on a count index | | `--csv` | off | Write the matrix as plain CSV instead of the default relaxed-PHYLIP format | | `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a `snp-*` value | | `--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 | Every value routes to one of two independent computations: ### Whole-index metrics | Value | Definition | |---|---| | `jaccard` | $D = 1 - \dfrac{\lvert A \cap B \rvert}{\lvert A \cup B \rvert}$ 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}}$ | `hamming` requires a presence/absence index; the others work on either index type. ### `snp-*` corrections Computed from the central-position SNP model (see "Central-position SNP model" below): a family is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base. These values require the sibling annex (`--sibling-annex`, below) and are, by default, computed exhaustively over every non-monomorphic family in the index; add `--subsample N` to bound the computation to approximately `N` families instead (see "Sampling at scale" below — the same flag `--pseudo-alignment`/`--sankoff` use, but optional here). For a genome pair, let $L$ be its total number of eligible loci (both genomes single-copy at that family), $p$ the raw proportion of substitutions among those loci, $P$/$Q$ the transition/transversion proportions, $Q_1$/$Q_2$ Kimura's two transversion categories (A↔C & G↔T vs. A↔T & C↔G), $P_1$/$P_2$ the purine (A↔G) / pyrimidine (C↔T) transition proportions, and $\pi_A,\pi_C,\pi_G,\pi_T$ the pair's pooled base frequencies. **`snp-raw`** $$d = p$$ **`snp-jc`** $$d = -\frac{3}{4}\ln\!\left(1-\frac{4p}{3}\right)$$ **`snp-k2p`** $$ \begin{aligned} a_1 &= 1-2P-Q \\ a_2 &= 1-2Q \\ d &= -\frac{1}{2}\ln a_1-\frac{1}{4}\ln a_2 \end{aligned} $$ **`snp-k81`** $$ \begin{aligned} a_1 &= 1-2P-2Q_1 \\ a_2 &= 1-2P-2Q_2 \\ a_3 &= 1-2Q_1-2Q_2 \\ d &= -\frac{1}{4}\left(\ln a_1+\ln a_2+\ln a_3\right) \end{aligned} $$ **`snp-f81`** $$ \begin{aligned} E &= 1-\left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right) \\ d &= -E\ln\!\left(1-\frac{p}{E}\right) \end{aligned} $$ **`snp-t92`** $$ \begin{aligned} g &= \pi_C+\pi_G \\ w &= 2g(1-g) \\ a_1 &= 1-\frac{P}{w}-Q \\ a_2 &= 1-2Q \\ d &= -w\ln a_1-\frac{1}{2}(1-w)\ln a_2 \end{aligned} $$ **`snp-tn93`** $$ \begin{aligned} g_R &= \pi_A+\pi_G \\ g_Y &= \pi_C+\pi_T \\ k_1 &= \frac{2\pi_A\pi_G}{g_R} \\ k_2 &= \frac{2\pi_C\pi_T}{g_Y} \\ k_3 &= 2\left(g_Rg_Y-\frac{\pi_A\pi_G\,g_Y}{g_R}-\frac{\pi_C\pi_T\,g_R}{g_Y}\right) \\ w_1 &= 1-\frac{P_1}{k_1}-\frac{Q}{2g_R} \\ w_2 &= 1-\frac{P_2}{k_2}-\frac{Q}{2g_Y} \\ w_3 &= 1-\frac{Q}{2g_Rg_Y} \\ d &= -k_1\ln w_1-k_2\ln w_2-k_3\ln w_3 \end{aligned} $$ **`snp-tv`** — transversions only, deliberately uncorrected: $$d = Q$$ `--gamma-shape ALPHA` applies to every value above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term in the formulas above is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (the same weight, same $x$). ### Output Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the first line, then one `labelvalue...` row per genome). With `--csv`, the format is instead a header row `genome,,,...` followed by one `