obikmer/UserDocMD/usage/phylo.md

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

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|auto	none	Rate-heterogeneity correction, for snp-* values that support it (see below). Either a fixed \alpha or auto/estimate to fit it from the data (see "Automatic \alpha estimation" 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).

Automatic \alpha estimation (--gamma-shape auto)

--gamma-shape auto (or the equivalent --gamma-shape estimate) fits \alpha from the index itself instead of requiring a user-supplied value, using a method-of-moments estimator computed once, from the same sampling pass that builds the pairwise substitution tally — no extra scan of the index.

The estimator pools substitution counts by partition rather than by genome pair: for partition i, let n_i be the total number of substitutions observed across every genome pair, and L_i the total number of eligible loci across every genome pair, in that partition. Define the partition's observed substitution rate:

R_i = \frac{n_i}{L_i}

Under a single shared substitution rate with no among-site heterogeneity, each R_i would vary only by Poisson sampling noise. Rate heterogeneity is modeled, as elsewhere in this correction, by a $\mathrm{Gamma}(\alpha,\alpha)$-distributed multiplicative rate (mean 1) shared by every locus in a partition — the classical Poisson–Gamma (negative-binomial) mixture. Under that model:

\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha}

where \mu is the pooled substitution rate across every partition. Weighting each partition's squared deviation by its own L_i removes the first (Poisson) term before attributing what's left to genuine rate heterogeneity:

\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad V = \frac{\sum_i L_i\,(R_i-\hat\mu)^2}{\sum_i L_i} \qquad \bar L = \frac{\sum_i L_i}{\text{number of partitions}} \hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L}

If the measured variance V doesn't exceed the Poisson floor \hat\mu/\bar L (no detectable over-dispersion across partitions — the data are consistent with a single shared rate), \alpha is left undefined: the correction is silently disabled for that run rather than applying a fabricated value, and a warning is logged. When an estimate is produced, it's logged at the info level before the distance matrix is computed.

Note: this is a method-of-moments estimator derived from the standard Poisson–Gamma relationship between substitution counts and gamma-distributed rate variation, applied per-partition — it is not part of Jin & Nei's (1990) original publication, which only defines the +Γ distance formula itself and, absent an estimate, recommends the fixed default \alpha = 1 (--gamma-shape 1) rather than proposing a way to estimate it from data. alpha < 1 indicates strong among-site rate heterogeneity (many near-invariant loci, a few fast ones); alpha growing large makes the correction converge to the uncorrected formula.

Output

Without -o, the matrix goes to stdout in relaxed-PHYLIP format (n on the first line, then one label<TAB>value... row per genome). With --csv, the format is instead a header row genome,<label1>,<label2>,... followed by one <label>,<value1>,<value2>,... row per genome, 6 decimals. Both formats are symmetric with a zero diagonal, except where noted below.

--exclude-genome, --min-shared-family

Option	Description
--exclude-genome LABEL	Exclude a genome (repeatable). Drops its row/column from the distance/shared-kmer matrix output, and removes it from the sampling used by --pseudo-alignment/--sankoff/a snp-* --distance value. Does not change the value computed for any remaining pair
--min-shared-family N	Auto-exclude, on top of --exclude-genome, any genome whose mean shared-family count against every other genome (see "Family Overlap" below) falls below N. Applies only to --pseudo-alignment/--sankoff/snp-* --distance — never to the whole-index metrics or their matrix/NJ/UPGMA output

Neighbor-Joining and UPGMA trees (--nj/--upgma) are always built from every genome in the index, regardless of --exclude-genome/--min-shared-family.

Central-position SNP model

Requires the sibling annex, built once per index:

Option	Description
--sibling-annex	Build (or rebuild) the sibling-count/minorant annex — prerequisite for every option in this section, and for a snp-* --distance value
--sibling-stats	Write <prefix>_siblings.csv: the family-size distribution, per genome and globally
--sibling-hist	Print the global family-size histogram (1-4 members) only
--family-overlap	Write <prefix>_family_overlap.csv: for every genome pair, how many variable families both genomes carry a call for
--shannon	Write <prefix>_entropy.csv: per-family Shannon entropy, one row per family, full unsampled scan
--pseudo-alignment	Write <prefix>_alignment.fasta: a SNP-only pseudo-alignment. Requires --subsample N

A family is eligible for a genome pair (i,j) only if both genomes carry exactly one of its observed forms (single-copy, unambiguous).

--sibling-stats

<prefix>_siblings.csv — family size = number of distinct central bases observed at a family (1-4).

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

Family Overlap

--family-overlap writes <prefix>_family_overlap.csv: header genome,<label1>,<label2>,..., one row per genome, cell [i][j] = number of variable families (family size ≥ 2) where both genome i and genome j carry a call. The diagonal is always 0. Every genome is written, unfiltered by --exclude-genome/--min-shared-family.

--min-shared-family N uses the mean of each genome's own row (excluding the diagonal) against this same matrix as its exclusion statistic. There is no universal value for N — inspect --family-overlap's own output to find where the real gap sits in a given genome collection before choosing a threshold.

--pseudo-alignment

<prefix>_alignment.fasta — one record per non-excluded genome, one column per variable family (family size ≥ 2). Each site is IUPAC-coded from the genome's presence mask at that family: a single observed form → the plain base; several forms → the matching IUPAC ambiguity code; no form → -.

Option	Default	Description
--subsample N	none (mandatory here)	Target number of families to sample
--free-loss	off	Treat a genome carrying none of a family's observed members as missing data (?) instead of -
--no-ambiguity	off	Treat a genome carrying more than one member of a family as missing data (?) instead of an IUPAC ambiguity code
--entropy MU	off (1.0 if only --entropy-sd is given)	Center of the entropy band to favor when sampling
--entropy-sd SIGMA	off (0.5 if only --entropy is given)	Width of that band

Sampling at scale: --subsample, --shannon, --entropy

--subsample, --free-loss, --no-ambiguity, --entropy/--entropy-sd are shared by --pseudo-alignment, --sankoff (and everything it implies: --tnt/--phyg/--iqtree), and a snp-* --distance value — one draw feeds all of them in a single invocation. --subsample is mandatory for --pseudo-alignment/--sankoff; for a snp-* --distance value it is optional (omitted means every non-monomorphic family in the index, not an approximation).

Combining --sankoff (or --tnt/--phyg/--iqtree) with a snp-* --distance value in the same command reuses that one draw for both — the distance and the Sankoff calibration/alignment are guaranteed to be computed from the identical set of sampled sites, never two independent samples, so the two outputs are directly comparable. This only holds within a single command; running them as two separate obikmer phylo invocations draws two independent samples even with the same flags.

Without --subsample, every variable family (family size ≥ 2) is used. With --subsample N, roughly N families are kept instead, drawn in proportion to how many candidate families each part of the index actually holds, so the sample stays representative of the whole index. If the index has fewer than N candidate families, --subsample has no effect.

--shannon: measuring how informative a family is

<prefix>_entropy.csv has one row per family visited:

Column	Meaning
layer	an internal index-layer identifier — stable within one run, not meaningful across indexes
family_idx	the family's position within that layer
entropy15	Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of {A,C,G,T}), genomes absent from the family excluded from the count
entropy4	Shannon entropy (bits) reduced to the 4 plain bases, kept alongside entropy15 for comparison
family_size	number of distinct central bases observed anywhere in the index for this family (2-4)
n_genomes_present	how many genomes the entropy was computed over

Run with --subsample N --shannon to get a bounded diagnostic sample instead of a full-index pass — useful for choosing --entropy/--entropy-sd values before a full run.

--entropy MU / --entropy-sd SIGMA: biasing the sample toward informative families

By default, --subsample draws families uniformly. With --entropy/--entropy-sd, each family's chance of being kept is instead weighted by how close its own entropy (entropy15) is to MU, using a Gaussian curve of width SIGMA — no hard cutoff. The filter activates as soon as either flag is given; the other defaults to 1.0/0.5. Combine with --subsample N (expect somewhat fewer than N families kept in practice) or use alone (a soft filter over the whole index, no size target).

The first phylo run on a given index that uses --entropy/--entropy-sd pays a one-time extra cost (every candidate family's entropy is computed once and saved alongside the index); later runs, even with different MU/SIGMA, reuse that saved data.

Sankoff calibration and phylogenetic exports

Option	Default	Description
--sankoff	off	Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment. Requires --subsample N
--sankoff-ratio-ceiling	0.5	Exclude genome pairs whose raw SNP ratio exceeds this value from the base-composition part of the calibration
--free-loss	off	Recode a family's non-detection as the ? missing-data symbol instead of an ordinary, costed state, throughout --sankoff 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)
--iqtree-min-freq	0.001	With --iqtree --free-loss: also treat as missing any state rarer than this in the alignment
--sankoff-cost-scale	100	Integer scaling factor applied to costs before rounding, for TNT/PhyG's integer-only cost commands

The 16-state model

Each family is a character with 16 possible states: one per subset of the 4 possible central bases (including the empty subset). Calibration combines a 5 \times 5 transition matrix over family cardinality (0-4 observed forms) and a 4 \times 4 base-substitution matrix from unambiguous single-copy loci, both restricted to genome pairs at or below --sankoff-ratio-ceiling, 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). It does not run any external tool.

With --free-loss, the empty state is recoded to ? (TNT/PhyG/IQ-TREE's own missing-data symbol) instead of an ordinary, costed 16th state — ? rather than -, since - still carries gap/indel semantics in these tools. --free-loss also zeroes the cardinality-transition cost between any two states, not just to/from the empty one: gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family.

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. 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 order the genomes appear in <prefix>_sankoff.fasta).

Output files

With -o/--output PREFIX, the relevant subset of the files below is written. Without -o, only the distance matrix is produced, on stdout. All matrices use genome labels as row/column headers, in index order.

Distance matrix

File	Written by	Format	Content
<prefix>_dist.phy	always, unless --csv	relaxed PHYLIP	the --distance matrix
<prefix>_dist.csv	--csv	CSV matrix	the --distance matrix, 6 decimals
<prefix>_shared.csv	--shared-kmers	CSV matrix	shared-kmer count per genome pair (integers)
<prefix>_nj.nwk	--nj	Newick	Neighbor-Joining tree
<prefix>_upgma.nwk	--upgma	Newick	UPGMA tree

CSV matrix layout (_dist.csv, _shared.csv, _family_overlap.csv): 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>_family_overlap.csv	--family-overlap	CSV matrix	variable families both genomes of a pair carry a call for
<prefix>_entropy.csv	--shannon	CSV table	per-family Shannon entropy, see "Sampling at scale" above
<prefix>_alignment.fasta	--pseudo-alignment	FASTA	SNP-only pseudo-alignment, IUPAC-coded

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, header carries an n_sites annotation
<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
<prefix>_iqtree_states.csv	--iqtree	CSV table	maps _iqtree.model/_iqtree.fasta's compact state symbols back to _sankoff_matrix.csv's alphabet

_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. One row per source state, one value per destination state, cost -\ln P(a,b), 4 decimals.

_sankoff_params.yaml:

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 — 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) — 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).

_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).

_sankoff.pg (--phyg) — script: read(prefasta:..., tcm:...) against _sankoff.fasta/_sankoff.tcm, a default 300s/4-instance search, report(...) writing <prefix>_sankoff.tre. Run with phyg from the output directory (the script uses relative file 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); with --iqtree-min-freq also set (the default), any state rarer than that threshold is folded into the same ? treatment, and non-informative columns are re-checked and dropped again. Run with:

iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO

_iqtree_states.csv (--iqtree) — one row per state actually kept in _iqtree.model/_iqtree.fasta (header iqtree_symbol,canonical_symbol,frequency): iqtree_symbol is the compact 0-9A-F symbol as written in those two files, canonical_symbol is the matching _sankoff_matrix.csv state, frequency is that state's empirical frequency at full precision. Under --free-loss, absent (0/?) is never a kept state, so it never appears here — nor does any state --iqtree-min-freq folded away for being too rare.

Rare states and --iqtree-min-freq

States that combine 3 or 4 central bases at once (IUPAC V/H/K.../N) are inherently rare, and can make iqtree3 itself numerically unstable ("Numerical underflow for lh-derivative" warnings). With --free-loss set, --iqtree-min-freq (default 0.001, one in a thousand) extends the missing-data treatment to any state below this frequency, not just absence. Check _iqtree_states.csv to see exactly which states survived and at what frequency; set --iqtree-min-freq 0 to keep every state that occurs at all. Has no effect without --free-loss.