From 52858dfa4aff1e63b5bf39971f3703c42c378441 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 12 Sep 2026 18:27:15 +0200 Subject: [PATCH] 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. --- _Sidebar.md | 2 +- ...phylogeny-kmer_based-kmer_set_distances.md | 71 +++++++++++++++++++ ...y-phylogeny-snp_based-central_snp_model.md | 2 +- usage-phylo.md | 6 +- 4 files changed, 76 insertions(+), 5 deletions(-) create mode 100644 theory-phylogeny-kmer_based-kmer_set_distances.md diff --git a/_Sidebar.md b/_Sidebar.md index 2ce9578..74b964f 100644 --- a/_Sidebar.md +++ b/_Sidebar.md @@ -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) diff --git a/theory-phylogeny-kmer_based-kmer_set_distances.md b/theory-phylogeny-kmer_based-kmer_set_distances.md new file mode 100644 index 0000000..0009be1 --- /dev/null +++ b/theory-phylogeny-kmer_based-kmer_set_distances.md @@ -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. diff --git a/theory-phylogeny-snp_based-central_snp_model.md b/theory-phylogeny-snp_based-central_snp_model.md index a737508..d78c605 100644 --- a/theory-phylogeny-snp_based-central_snp_model.md +++ b/theory-phylogeny-snp_based-central_snp_model.md @@ -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. diff --git a/usage-phylo.md b/usage-phylo.md index 09b815d..80ad997 100644 --- a/usage-phylo.md +++ b/usage-phylo.md @@ -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`.