Refactor phylogenetic models and update usage documentation

This commit removes entries related to Kmer-based and SNP-based phylogenetic methods, refining the definitions for SNP models, and introducing new documentation for k-mer based set-based distances. Additionally, new configuration options and exclusion flags are added to clarify tree generation behavior.
Eric Coissac committed 2026-09-12 18:27:37 +02:00
1 parent 19f4ec2567
commit 52858dfa4a
4 files changed
+76 -5

No files matched your search

+1 -1
@@ -12,7 +12,7 @@
- [Low-complexity kmer filter](theory-kmer_indexing-entropy_filter)
### Phylogeny
#### Kmer-based
- [Distance metrics](theory-phylogeny-kmer_based-distance_metrics)
- [Whole k-mer set distances](theory-phylogeny-kmer_based-kmer_set_distances)
#### SNP-based
- [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model)
- [Gamma rate heterogeneity](theory-phylogeny-snp_based-gamma_rate_heterogeneity)
+71
@@ -0,0 +1,71 @@
# Whole k-mer set distances
The k-mer composition of a genome provides a representation of its sequence content that can be compared directly between genomes. Distances in this section measure the dissimilarity between two genomes from their complete collection of k-mers.
Depending on the information retained about the k-mers, two complementary types of comparison can be considered. **Set-based distances** use only the presence or absence of each distinct k-mer, whereas **count-based distances** also account for how frequently each k-mer occurs. Count-based distances can in turn be expressed either from raw k-mer counts or from the corresponding relative frequencies.
Let $A$ and $B$ denote the k-mer sets of two genomes. For each distinct k-mer $i$, let $c_i^A$ and $c_i^B$ denote its counts in the two genomes, and let
$$p_i^A=\frac{c_i^A}{\sum_j c_j^A}, \qquad p_i^B=\frac{c_i^B}{\sum_j c_j^B}$$
denote the corresponding relative frequencies.
These distances can then be classified along two independent axes: **the information retained about k-mer occurrences** (presence/absence versus abundance) and **whether the resulting dissimilarity is a metric**, i.e. whether it satisfies the triangle inequality.
## Presence/absence-based distances
When only presence or absence is retained, a genome is reduced to the set of k-mers it contains, and a distance can only compare which k-mers the two genomes share, not how often each occurs.
Two distances built this way are true metrics, meaning they satisfy the triangle inequality — the distance from genome $A$ to genome $C$ can never exceed the sum of the distances from $A$ to $B$ and from $B$ to $C$, whichever genome $B$ is chosen as an intermediate. The **Jaccard distance** compares the two k-mer sets through the size of their intersection relative to their union:
$$D = 1 - \frac{\lvert A \cap B \rvert}{\lvert A \cup B \rvert}$$
The **Hamming distance** instead counts, k-mer by k-mer over the whole k-mer space, how many positions disagree on presence or absence between the two genomes:
$$D = \sum_i \mathbb{1}[a_i \ne b_i]$$
without normalizing by the number of positions compared. Because it is defined directly on presence/absence, Hamming distance cannot be computed from counts the way the other measures below can.
A third distance, **Mash**, is derived from the Jaccard distance rather than computed independently. It approximates a per-site mutation rate between the two genomes from their Jaccard distance $J$ and the k-mer size $k$:
$$D = -\frac{1}{k} \ln\!\left(\frac{2(1-J)}{2-J}\right)$$
(clamped to 1.0 when the Jaccard distance reaches its maximum). This transformation is monotonic — it never reorders which of two genome pairs is closer or farther — but a monotonic transformation of a metric is not guaranteed to remain a metric itself, and whether the Mash distance satisfies the triangle inequality has not been established here. It should be treated as a useful evolutionary approximation rather than a proven metric.
## Abundance-based distances
When k-mer counts are retained, either as raw counts $c_i$ or as the relative frequencies $p_i$ introduced above, distances can additionally reflect how much more or less abundant a shared k-mer is between the two genomes, not merely whether it is present in both.
The **Euclidean distance**, computed on either raw counts or relative frequencies, is the ordinary straight-line distance between the two genomes' count vectors:
$$D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$$
Euclidean distance is always a metric, on counts as on frequencies (`relfreq-euclidean`).
The **Hellinger distance** is built the same way, but on the square roots of the relative frequencies rather than the frequencies themselves:
$$D = \frac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$$
and stays bounded between 0 and 1. Because it is, by construction, a rescaled Euclidean distance between the $\sqrt{p}$ vectors, it is itself a metric — as is `hellinger-euclidean`, its unnormalized variant $D = \sqrt{2} \times D_{\text{hellinger}}$, since multiplying a metric by a positive constant cannot break the triangle inequality.
The **Bray-Curtis dissimilarity** compares two count vectors through how much of their combined total is *not* shared:
$$D = 1 - \frac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$$
Unlike the distances above, Bray-Curtis computed on raw counts is *not* a metric in general — it can violate the triangle inequality, a well-documented property of this measure rather than an approximation error. The reason becomes clear once the formula is rewritten, using $a+b-2\min(a,b) = \lvert a-b \rvert$, as
$$D = \frac{\sum_i \lvert c_i^A - c_i^B \rvert}{\sum_i c_i^A + \sum_i c_i^B}$$
The denominator is the combined total abundance of both genomes, and it changes from one genome pair to the next whenever their total k-mer counts differ — it is precisely this shifting denominator that breaks the triangle inequality.
That reasoning also shows exactly when Bray-Curtis *does* behave as a metric: whenever every count vector being compared shares the same fixed total $S$, the denominator becomes the constant $2S$ for every pair, and the formula reduces to
$$D = \frac{1}{2S}\sum_i \lvert c_i^A - c_i^B \rvert$$
a positive multiple of the $L^1$ distance, which is always a metric. Relative frequencies are exactly such a case: every genome's relative frequencies sum to $S = 1$ by construction, so **`relfreq-bray-curtis`** reduces to
$$D = \frac{1}{2}\sum_i \lvert p_i^A - p_i^B \rvert$$
the standard total variation distance between two probability distributions, a well-known metric. `bray-curtis` on raw counts and `relfreq-bray-curtis` on relative frequencies are therefore the same formula applied to two different kinds of vectors, and only one of the two — the one where every vector has the same total — is guaranteed to behave as a proper distance.
See [`phylo`](usage-phylo) for how to select a distance and the resulting output formats.
+1 -1
@@ -2,7 +2,7 @@
## Family
A **family** is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base — the odd kmer length guarantees a single, well-defined central position (see [Kmers](theory-kmer_indexing-kmers)). A family is **eligible** for a genome pair $(i,j)$ only if both genomes carry exactly one of its observed forms (single-copy, unambiguous) — this excludes multi-copy and absent loci from the comparison rather than folding them into an undifferentiated "not identical" bucket the way a whole-index metric would.
A **family** is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base — the odd kmer length guarantees a single, well-defined central position (see [Kmers](theory-kmer_indexing-kmers)). A family is **eligible** for a genome pair $(i,j)$ only if both genomes carry exactly one of its observed forms (single-copy, unambiguous) — this excludes multi-copy and absent loci from the comparison rather than folding them into an undifferentiated "not identical" bucket the way a whole-index distance would.
Conditioning on eligibility this way restricts every comparison to loci that are directly, positively confirmed comparable in both genomes: a locus never enters the statistic because of a genuinely absent homologous region, a diverged paralogous copy, or a genome-size asymmetry — only because a real single-copy substitution (or lack of one) was observed at flanks confirmed intact in both genomes.
+3 -3
@@ -20,12 +20,12 @@ obikmer phylo INDEX [OPTIONS]
| `--gamma-shape ALPHA\|auto` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see [Gamma rate heterogeneity](theory-phylogeny-snp_based-gamma_rate_heterogeneity)). Either a fixed $\alpha$ or `auto`/`estimate` to fit it from the data. 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 |
| `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index distance, 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** (`jaccard`, `mash`, `hamming`, `bray-curtis`, `relfreq-bray-curtis`, `euclidean`, `relfreq-euclidean`, `hellinger`, `hellinger-euclidean` — see [Distance metrics](theory-phylogeny-kmer_based-distance_metrics) for the definitions; `hamming` requires a presence/absence index, the others work on either index type), or **`snp-*` corrections** (`snp-raw`, `snp-jc`, `snp-k2p`, `snp-k81`, `snp-f81`, `snp-t92`, `snp-tn93`, `snp-tv` — see the [central-position SNP model](theory-phylogeny-snp_based-central_snp_model) for the definitions).
Every value routes to one of two independent computations: **whole-index distances** (`jaccard`, `mash`, `hamming`, `bray-curtis`, `relfreq-bray-curtis`, `euclidean`, `relfreq-euclidean`, `hellinger`, `hellinger-euclidean` — see [Whole k-mer set distances](theory-phylogeny-kmer_based-kmer_set_distances) for the definitions; `hamming` requires a presence/absence index, the others work on either index type), or **`snp-*` corrections** (`snp-raw`, `snp-jc`, `snp-k2p`, `snp-k81`, `snp-f81`, `snp-t92`, `snp-tn93`, `snp-tv` — see the [central-position SNP model](theory-phylogeny-snp_based-central_snp_model) for the definitions).
`snp-*` 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).
@@ -38,7 +38,7 @@ Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the fir
| 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 |
| `--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 distances 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`.