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

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.

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

On a large index (billions of families), building a full pseudo-alignment or fully calibrating the Sankoff model is not just slow — it may not fit in the time you have. --subsample bounds the work to a fixed number of families; --shannon reports how informative each family is; --entropy/--entropy-sd bias which families get kept toward the informative ones instead of choosing uniformly at random.

Option	Default	Description
--subsample N	none (keep everything)	Cap the number of variable families (family size ≥ 2) retained, to approximately N
--shannon	off	Write <prefix>_shannon.csv: per-family Shannon entropy, one row per family
--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

--subsample/--entropy/--entropy-sd affect every option that scans variable families: --snp, --family-overlap, --shannon, and the whole Sankoff pipeline (--sankoff/--tnt/--phyg/--iqtree, next section) — all of them draw from the same selection of families in one invocation, so the Sankoff calibration and the pseudo-alignment it calibrates always describe the same sites, and --family-overlap's counts stay consistent with --snp's columns. --raw-snp-distance/--raw-snp-counts are not affected — they always scan every family, since their p-distance estimate is a whole-index statistic, not something that benefits from being restricted to a sample.

--subsample N

Without --subsample, every variable family (family size ≥ 2, i.e. every family where at least one genome differs from the rest) is used. With --subsample N, roughly N families are kept instead, chosen at random but in proportion to how many candidate families each part of the index actually holds — so the sample stays representative of the whole index, not skewed toward whichever part happens to be scanned first. If the index has fewer than N candidate families in the first place, --subsample has no effect: everything is kept.

--subsample trades completeness for speed: --snp's alignment gets fewer columns, --sankoff's calibration is based on fewer observations, but the resolution work (the expensive part of a phylo run on a large index) scales with N instead of with the index's true size. Pick N as large as your time budget allows — a few hundred thousand to a few million families is usually enough for the transition-probability estimates in --sankoff's calibration to stabilize; a smaller N speeds up exploratory runs.

--shannon: measuring how informative a family is

Not every variable family is equally useful for a tree: a family that differs in only one genome out of a thousand carries very little signal, and one where the pattern looks essentially random across genomes may be too saturated (multiple substitutions have overwritten the original signal) to carry real information either. --shannon quantifies this with the Shannon entropy (in bits) of each family's states across the genomes that carry it — low entropy means "almost everyone agrees" (an invariant or near-invariant family, phylogenetically shallow), while entropy near the ceiling for a 4-state character means "close to a random draw between the possible bases" (saturated).

<prefix>_shannon.csv (or shannon.csv without -o) 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} — the same alphabet --sankoff's 16-state model uses), genomes absent from the family excluded from the count
entropy4	Shannon entropy (bits) reduced to the 4 plain bases, kept alongside entropy15 for comparison — a genome carrying more than one base at once counts once per base, so this can differ from entropy15
family_size	number of distinct central bases observed anywhere in the index for this family (2–4, since monomorphic families aren't visited)
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 to inspect the entropy distribution and decide reasonable --entropy/--entropy-sd values (see below) before committing to a full run.

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

By default, --subsample draws families uniformly — every candidate family has the same chance of being kept, regardless of how informative it actually is. --entropy/--entropy-sd change that: instead of a uniform draw, each family's chance of being kept is weighted by how close its own entropy (the entropy15 value --shannon reports) is to MU, using a bell-shaped (Gaussian) curve of width SIGMA — a family with entropy exactly MU is the most likely to be kept, and the chance falls off smoothly the further its entropy is from MU, with no hard cutoff (a few families outside the target band can still get in, just less often).

The filter activates as soon as either --entropy or --entropy-sd is given; whichever one you don't set defaults to 1.0/0.5. It can be combined with --subsample N (the target count is still approximately N, now biased toward the entropy band instead of uniform — expect somewhat fewer than N families in practice, since low-weight families are dropped rather than replaced) or used alone (--entropy without --subsample: a soft entropy 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 has to be computed once, up front, and is then saved alongside the index). Every following run — even with different MU/SIGMA values — reuses that saved data and stays fast.

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 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
<prefix>_shannon.csv	--shannon	CSV table	per-family Shannon entropy, see "Sampling at scale" above

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