Refactor Kmer indexing and add phylogenetic analysis features
This commit introduces a complete overhaul of the Kmer indexing architecture, including new support for 2-bit encoding, canonical super-kmer handling, and a partitioned processing pipeline. Additionally, it adds comprehensive phylogenetic capabilities, including multiple distance metrics, SNP correction models, and models for rate heterogeneity and sampling design.
Eric Coissac committed
1 parent
0b79a8ee74
commit
97fe82159a
16 files changed
+423
-98
No files matched your search
+4
-4
@@ -6,11 +6,11 @@ All functionality is exposed through a single binary, `obikmer`, organized as su
|
||||
|
||||
## Core principles
|
||||
|
||||
- Kmers are of fixed, odd length $k$, chosen at index-construction time in the range $[11, 31]$ (see [Kmers and super-kmers](theory-kmers_and_superkmers)).
|
||||
- Each kmer fits in a 64-bit word using a 2-bit-per-base encoding (see [DNA encoding](theory-encoding)).
|
||||
- Kmers are of fixed, odd length $k$, chosen at index-construction time in the range $[11, 31]$ (see [Kmers](theory-kmer_indexing-kmers)).
|
||||
- Each kmer fits in a 64-bit word using a 2-bit-per-base encoding (see [DNA encoding](theory-kmer_indexing-encoding)).
|
||||
- Kmers are handled in **canonical form** ($\text{canonical}(kmer) = \min(kmer, \text{revcomp}(kmer))$), making counting strand-independent.
|
||||
- Sequences are decomposed into **super-kmers** before storage, anchored on a hash-selected **minimizer** (see [Minimizer selection](theory-minimizer_selection)), then routed to one of several **partitions** for parallel, memory-bounded processing (see [Partitioning and indexing architecture](theory-indexing_architecture)).
|
||||
- Low-complexity kmers can be filtered out at index-construction time using an entropy-based score (see [Low-complexity kmer filter](theory-entropy_filter)).
|
||||
- Sequences are decomposed into **super-kmers** before storage, anchored on a hash-selected **minimizer** (see [Minimizer selection](theory-kmer_indexing-minimizer_selection)), then routed to one of several **partitions** for parallel, memory-bounded processing (see [Partitioning and indexing architecture](theory-kmer_indexing-indexing_architecture)).
|
||||
- Low-complexity kmers can be filtered out at index-construction time using an entropy-based score (see [Low-complexity kmer filter](theory-kmer_indexing-entropy_filter)).
|
||||
|
||||
## Commands
|
||||
|
||||
|
||||
+15
-5
@@ -3,11 +3,21 @@
|
||||
- [Home](Home)
|
||||
- [Installation](installation)
|
||||
## Theory
|
||||
- [Kmers and super-kmers](theory-kmers_and_superkmers)
|
||||
- [DNA encoding](theory-encoding)
|
||||
- [Low-complexity kmer filter](theory-entropy_filter)
|
||||
- [Minimizer selection](theory-minimizer_selection)
|
||||
- [Partitioning and indexing architecture](theory-indexing_architecture)
|
||||
### Kmer indexing
|
||||
- [DNA encoding](theory-kmer_indexing-encoding)
|
||||
- [Kmers](theory-kmer_indexing-kmers)
|
||||
- [Minimizer selection](theory-kmer_indexing-minimizer_selection)
|
||||
- [Super-kmers](theory-kmer_indexing-superkmers)
|
||||
- [Partitioning and indexing architecture](theory-kmer_indexing-indexing_architecture)
|
||||
- [Low-complexity kmer filter](theory-kmer_indexing-entropy_filter)
|
||||
### Phylogeny
|
||||
#### Kmer-based
|
||||
- [Distance metrics](theory-phylogeny-kmer_based-distance_metrics)
|
||||
#### SNP-based
|
||||
- [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model)
|
||||
- [Gamma rate heterogeneity](theory-phylogeny-snp_based-gamma_rate_heterogeneity)
|
||||
- [Sankoff model](theory-phylogeny-snp_based-sankoff_model)
|
||||
- [Sampling design](theory-phylogeny-snp_based-sampling_design)
|
||||
## Usage
|
||||
- [superkmer](usage-superkmer)
|
||||
- [index](usage-index_command)
|
||||
|
||||
+2
-2
@@ -1,6 +1,6 @@
|
||||
# Architecture notes for advanced use
|
||||
|
||||
This page describes execution-level behavior relevant to sizing and running `obikmer` on large datasets or multi-socket machines. It complements the [index format](formats-index_layout) and [theory](theory-indexing_architecture) pages.
|
||||
This page describes execution-level behavior relevant to sizing and running `obikmer` on large datasets or multi-socket machines. It complements the [index format](formats-index_layout) and [theory](theory-kmer_indexing-indexing_architecture) pages.
|
||||
|
||||
## Sequence invariant
|
||||
|
||||
@@ -8,7 +8,7 @@ Every input sequence is treated purely as a compact representation of a set of o
|
||||
|
||||
- Only the `A`/`C`/`G`/`T` alphabet (case-insensitive) is recognized; a sequence is cut at any other character (including IUPAC ambiguity codes), so runs containing them are not represented in the index.
|
||||
- Sequences are internally processed in chunks of at most 256 nucleotides; a chunk shorter than k is dropped. This is invisible to the user beyond the ACGT-only, minimum-length-k constraints above.
|
||||
- Kmers are always handled in canonical form (see [DNA encoding](theory-encoding)), so the tool is strand-agnostic throughout: a kmer and its reverse complement are always the same entry.
|
||||
- Kmers are always handled in canonical form (see [DNA encoding](theory-kmer_indexing-encoding)), so the tool is strand-agnostic throughout: a kmer and its reverse complement are always the same entry.
|
||||
|
||||
## Index dimensioning
|
||||
|
||||
|
||||
+2
-2
@@ -2,9 +2,9 @@
|
||||
|
||||
## Construction pipeline
|
||||
|
||||
Building an index ([`index`](usage-index_command)) proceeds through a fixed sequence of phases, each operating independently per partition (see [Partitioning and indexing architecture](theory-indexing_architecture)):
|
||||
Building an index ([`index`](usage-index_command)) proceeds through a fixed sequence of phases, each operating independently per partition (see [Partitioning and indexing architecture](theory-kmer_indexing-indexing_architecture)):
|
||||
|
||||
1. **Scatter.** A single streaming pass over the input. Each sequence fragment is cut at non-ACGT bases, passed through the low-complexity entropy filter (see [Low-complexity kmer filter](theory-entropy_filter)), and any resulting segment shorter than k is dropped. Surviving segments are decomposed into super-kmers, canonicalized, and routed by `hash(minimizer) mod n_partitions` into one file per partition.
|
||||
1. **Scatter.** A single streaming pass over the input. Each sequence fragment is cut at non-ACGT bases, passed through the low-complexity entropy filter (see [Low-complexity kmer filter](theory-kmer_indexing-entropy_filter)), and any resulting segment shorter than k is dropped. Surviving segments are decomposed into super-kmers, canonicalized, and routed by `hash(minimizer) mod n_partitions` into one file per partition.
|
||||
2. **Dereplication.** Within each partition, identical super-kmer sequences are merged and their occurrence counts summed. This count is per super-kmer, not per kmer — a kmer's true abundance is the sum of the counts of every super-kmer containing it.
|
||||
3. **Exact counting.** Every kmer in every dereplicated super-kmer is enumerated and its exact total count computed. A per-genome kmer frequency spectrum is produced at this stage.
|
||||
4. **Quorum filtering.** Kmers outside the `--min-abundance`/`--max-abundance` range are dropped, and super-kmers are recompacted around the surviving kmer set.
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
# DNA encoding
|
||||
|
||||
## 2-bit nucleotide encoding
|
||||
|
||||
Every nucleotide is encoded on 2 bits, most-significant-bit first within each word:
|
||||
|
||||
| Base | Encoding |
|
||||
| --- | --- |
|
||||
| A | `00` |
|
||||
| C | `01` |
|
||||
| G | `10` |
|
||||
| T | `11` |
|
||||
|
||||
The Watson-Crick complement of a base is its bitwise NOT on 2 bits: $\text{complement}(base) = \lnot base \mathbin{\&} \texttt{0b11}$.
|
||||
|
||||
## Kmer encoding
|
||||
|
||||
A kmer of length $k$ ($k \le 31$) fits in a single 64-bit word. The first nucleotide occupies the two most significant bits, each following nucleotide occupies the next two bits, and unused low-order bits are zero. Extracting nucleotide i (0-indexed from the 5′ end) is a shift-and-mask operation.
|
||||
|
||||
Reverse complement is computed by bit manipulation directly on the packed word, without any lookup table: complement every base, reverse the byte order, then reverse the order of 2-bit groups within each byte in two more passes, and finally realign the result to the most-significant bits.
|
||||
|
||||
## Canonical form
|
||||
|
||||
The canonical form of a packed 2-bit sequence is the lexicographic minimum of the sequence and its reverse complement:
|
||||
|
||||
$$\text{canonical}(x) = \min\big(x,\ \text{revcomp}(x)\big)$$
|
||||
|
||||
This is a single operation on the packed representation, so it applies identically regardless of what `x` represents — a kmer (see [Kmers](theory-kmer_indexing-kmers)), an m-mer (see [Minimizer selection](theory-kmer_indexing-minimizer_selection)), or a super-kmer (see [Super-kmers](theory-kmer_indexing-superkmers)). Using the canonical form halves the relevant space and makes counting/selection strand-independent: a sequence and its reverse complement are always treated as the same entity, regardless of which DNA strand was sequenced.
|
||||
@@ -0,0 +1,36 @@
|
||||
# Low-complexity kmer filter
|
||||
|
||||
Low-complexity kmers (homopolymer runs, tandem repeats) can dominate an index without carrying useful information. `obikmer` detects and excludes them during index construction using a normalized Shannon entropy score.
|
||||
|
||||
## Sub-word frequencies
|
||||
|
||||
For a kmer of length $k$ and a sub-word size $ws$ ($1 \le ws \le ws_{\max}$, default $ws_{\max} = 6$), the kmer is decomposed into its $k - ws + 1$ overlapping sub-words of length $ws$ by sliding a window across it. Each sub-word is tallied under its raw 2-bit-packed value, with no canonicalization.
|
||||
|
||||
## Corrected Shannon entropy
|
||||
|
||||
Let $f_j$ be the observed count of raw sub-word $j$, and $n_{\text{words}} = k - ws + 1$ the total number of sub-words. The entropy is:
|
||||
|
||||
$$H_{\text{corr}} = \log(n_{\text{words}}) - \frac{1}{n_{\text{words}}} \sum_j f_j \log f_j$$
|
||||
|
||||
## Small-sample correction
|
||||
|
||||
Because only $n_{\text{words}}$ sub-words are observed among up to $4^{ws}$ possible values, the achievable maximum entropy $H_{\max}$ is bounded below $\log(4^{ws})$ for small samples. $H_{\max}$ is computed from the most uniform integer distribution achievable with $n_{\text{words}}$ observations over $4^{ws}$ categories. The normalized entropy is:
|
||||
|
||||
$$\hat{H}(ws) = \frac{H_{\text{corr}}}{H_{\max}} \in [0, 1]$$
|
||||
|
||||
A value near 0 indicates low complexity (e.g. a homopolymer run); near 1 indicates high complexity, characteristic of a random sequence.
|
||||
|
||||
## Final score
|
||||
|
||||
The filter evaluates $\hat{H}(ws)$ for every word size from 1 to ws\_max and keeps the minimum:
|
||||
|
||||
$$\text{entropy}(kmer) = \min_{ws=1}^{ws_{\max}} \hat{H}(ws)$$
|
||||
|
||||
Taking the minimum across word sizes ensures that repetition at any scale is detected: a homopolymer is caught at $ws=1$, a dinucleotide repeat at $ws=2$, and so on. A kmer is rejected if its entropy score falls below a threshold $\theta$ (default 0.7), a configurable collection parameter.
|
||||
|
||||
## Properties
|
||||
|
||||
The entropy score depends only on the kmer sequence itself, not on where or how many times it occurs:
|
||||
|
||||
- **Orientation invariance**: a kmer and its reverse complement always receive the same score.
|
||||
- **Context independence**: a given kmer is always accepted or always rejected, regardless of which genome or read it appears in. The filter defines a fixed partition of the kmer space into low-complexity and valid kmers.
|
||||
@@ -0,0 +1,26 @@
|
||||
# Partitioning and indexing architecture
|
||||
|
||||
An index is split into a fixed number of **partitions**, each handling an independent, disjoint slice of the kmer space. Partitioning keeps the working set of each stage small enough to process efficiently and enables parallel construction and querying.
|
||||
|
||||
## Routing
|
||||
|
||||
The canonical minimizer of a super-kmer (see [Minimizer selection](theory-kmer_indexing-minimizer_selection)) is hashed to produce a $p$-bit routing value that selects the destination partition — the low $p$ bits of the minimizer's hash:
|
||||
|
||||
canonical minimizer → hash(minimizer) → p-bit value → partition index
|
||||
|
||||
$$\text{partition} = H(\text{minimizer}) \bmod 2^p$$
|
||||
|
||||
Within a partition, kmers are indexed as plain values via a minimal perfect hash function (see [On-disk storage](formats-index_layout)); the minimizer plays no further role once a super-kmer has reached its partition.
|
||||
|
||||
## Parameter guidance
|
||||
|
||||
Even though $H$ already makes minimizer values well-distributed (see [Minimizer selection](theory-kmer_indexing-minimizer_selection)), choosing $p$ well below the number of bits available in the minimizer ($2m$) leaves a comfortable entropy margin, provided the number of distinct minimizers actually observed is much larger than the number of partitions.
|
||||
|
||||
| Minimizer size $m$ | Minimizer bits ($2m$) | Typical partition-index bits $p$ | Partitions |
|
||||
| --- | --- | --- | --- |
|
||||
| 9 | 18 | 6–8 | 64–256 |
|
||||
| 11 | 22 | 8–10 | 256–1 024 |
|
||||
| 13 | 26 | 10–12 | 1 024–4 096 |
|
||||
| 15 | 30 | 10–14 | 1 024–16 384 |
|
||||
|
||||
The number of partitions must satisfy $p \le 2m$, and in practice $p$ is chosen well below that bound to leave a comfortable entropy margin. For $k=31$, $m=13$, $p=10$ (1024 partitions), partition load is well balanced on real genomic data.
|
||||
@@ -0,0 +1,6 @@
|
||||
# Kmers
|
||||
|
||||
A **kmer** is a DNA subsequence of fixed length $k$. Two constraints apply to $k$, both enforced when a command starts (an invalid value exits immediately with an error):
|
||||
|
||||
- $k \in [11, 31]$: long enough to be specific, short enough to fit in a single 64-bit word at 2 bits/base ($k \le 32$ is the hard limit; $k < 11$ gives insufficient specificity).
|
||||
- $k$ **is odd**: an odd-length sequence can never equal its own reverse complement, so the two orientations of any kmer are always distinct. This is required for the canonical form (see [DNA encoding](theory-kmer_indexing-encoding)) to be well defined.
|
||||
@@ -0,0 +1,125 @@
|
||||
# Minimizer selection
|
||||
|
||||
## Definition
|
||||
|
||||
The **minimizer** of a k-mer is the canonical form of the m-mer that is smallest according to a chosen ordering among the $k-m+1$ overlapping m-mers contained in the k-mer, with $m<k$ ([Roberts et al. 2004](#ref-Roberts2004-rz)). Canonical form is the same operation defined for a full kmer (see [DNA encoding](theory-kmer_indexing-encoding)), applied here to an m-mer $s$ instead: $s^c = \min_{\mathrm{lex}}\left(s,\operatorname{RC}(s)\right)$.
|
||||
|
||||
In **OBIkmer**, canonical m-mers are ordered according to a deterministic hash function $H$ applied to their canonical form:
|
||||
|
||||
$$s_1^c <_h s_2^c \quad\Longleftrightarrow\quad H(s_1^c) < H(s_2^c).$$
|
||||
|
||||
Thus, for a k-mer $K$, let
|
||||
|
||||
$$M(K)=(s_1,\ldots,s_{k-m+1})$$
|
||||
|
||||
be its sequence of overlapping m-mers. The minimizer is
|
||||
|
||||
$$\operatorname{minimizer}(K) = s_j^c, \qquad j=\underset{i}{\operatorname{argmin}}\;H(s_i^c).$$
|
||||
|
||||
In other words, the hash function defines the ordering of canonical m-mers, while the minimizer itself is the canonical form of the m-mer selected by that ordering.
|
||||
|
||||
Applying the minimizer to consecutive k-mers partitions a sequence into super-kmers (see [Super-kmers](theory-kmer_indexing-superkmers)).
|
||||
|
||||
## Hash-based minimizer ordering
|
||||
|
||||
`obikmer` selects minimizers by hash order rather than plain lexicographic order. Ordering m-mers lexicographically on their 2-bit encoding systematically favors AT-rich m-mers (an all-A m-mer always encodes to 0), which causes low-complexity regions to dominate as minimizers and produces unbalanced partitions ([Golan & Shur 2025](#ref-Golan2025-xf); [Kille et al. 2023](#ref-Kille2023-px); [Pan & Reinert 2024](#ref-Pan2024-hb); [Zheng et al. 2020](#ref-Zheng2020-ji); [2021](#ref-Zheng2021-cc)).
|
||||
|
||||
Instead, a well-distributed hash function $H$ is applied to the canonical (lexicographically minimal) form of each m-mer, and the m-mer with the smallest $H$ value wins. Because $H$ is a bijective mixing function with good avalanche properties ([Steele et al. 2014](#ref-Steele2014-cz)), its output behaves approximately as a composition-independent pseudorandom ordering of distinct m-mers. Under the usual random-hash assumption, each distinct m-mer in a window has approximately the same probability of being the minimum, independently of nucleotide composition.
|
||||
|
||||
The canonical form used as input to $H$ is still the lexicographic minimum of forward/reverse-complement — hashing is applied on top of it, not used to redefine it. Defining canonicity by hash value instead would bias the *distribution of hash values themselves* toward small values (the minimum of two independent hashes is not uniformly distributed), reintroducing a bias one layer down.
|
||||
|
||||
### Hash function
|
||||
|
||||
The hash function (H) applies Stafford's 64-bit mixing function, variant 13 ([Stafford 2011](#ref-Stafford2011)), to the m-mer encoding XORed with a seed $\sigma$:
|
||||
|
||||
```text
|
||||
H(x):
|
||||
x ← x ⊕ σ
|
||||
x ← x ⊕ (x >> 30)
|
||||
x ← x × 0xbf58476d1ce4e5b9
|
||||
x ← x ⊕ (x >> 27)
|
||||
x ← x × 0x94d049bb133111eb
|
||||
return x ⊕ (x >> 31)
|
||||
```
|
||||
|
||||
The hash uses a fixed non-zero seed
|
||||
|
||||
$$ \sigma=\texttt{0x9e3779b97f4a7c15}. $$
|
||||
|
||||
The choice of $\sigma$ is not arbitrary. Low-complexity m-mers (homopolymers and short tandem repeats) are disproportionately abundant in real genomes. If one of them happened to be the global argmin of $H$, it could therefore occur as the minimizer in far more windows than expected from the composition-uniform behavior described above. This is not a bias of the mixing function itself, but a consequence of the repeated occurrence of the same input in real sequence data.
|
||||
|
||||
In particular, with $\sigma=0$, the all-A m-mer is encoded as $0$, and
|
||||
|
||||
$$ \operatorname{Mix13}(0)=0, $$
|
||||
|
||||
This makes it the global argmin of the hash function. The non-zero seed eliminates the particular pathological correspondence between the all-A m-mer and the zero output of the mixer.
|
||||
|
||||
Exhaustive enumeration of all $4^m$ m-mers confirms that, with $\sigma=\texttt{0x9e3779b97f4a7c15}$, the global argmin is neither a homopolymer nor a periodic repeat for any tested (m):
|
||||
|
||||
| $m$ | argmin (canonical) | decoded sequence | minimal period |
|
||||
| --- | --- | --- | --- |
|
||||
| 3 | 16 | `CAA` | 3 |
|
||||
| 5 | 78 | `ACATG` | 5 |
|
||||
| 7 | 5512 | `CCCGAGA` | 7 |
|
||||
| 9 | 108760 | `CGGGATCGA` | 9 |
|
||||
| 11 | 179014 | `AAGGTGTCACG` | 11 |
|
||||
| 13 | 33759044 | `GAAATACTTCACA` | 13 |
|
||||
| 15 | 29869313 | `AACTACTTACCAAAC` | 15 |
|
||||
|
||||
When the minimal period of the sequence is $m$, as observed in the table above, the sequence is primitive: it cannot be represented as repetitions of a shorter sequence.
|
||||
|
||||
See [Partitioning and indexing architecture](theory-kmer_indexing-indexing_architecture) for how this hash routes a super-kmer to a partition.
|
||||
|
||||
## Bibliography
|
||||
|
||||
<div id="refs" class="references csl-bib-body hanging-indent">
|
||||
|
||||
<div id="ref-Golan2025-xf" class="csl-entry">
|
||||
|
||||
Golan, S. & Shur, A.M. (2025). <a href="https://doi.org/10.1007/978-3-031-82670-2\_25">Expected density of random minimizers</a>. In: <span style="font-style: italic;">Lecture Notes in Computer Science</span>, Lecture Notes in Computer Science. Springer Nature Switzerland, Cham, pp. 347–360.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Kille2023-px" class="csl-entry">
|
||||
|
||||
Kille, B., Garrison, E., Treangen, T.J. & Phillippy, A.M. (2023). <a href="https://doi.org/10.1093/bioinformatics/btad512">Minmers are a generalization of minimizers that enable unbiased local Jaccard estimation</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 39.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Pan2024-hb" class="csl-entry">
|
||||
|
||||
Pan, C. & Reinert, K. (2024). <a href="https://doi.org/10.1093/bioinformatics/btae045">A simple refined DNA minimizer operator enables 2-fold faster computation</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 40.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Roberts2004-rz" class="csl-entry">
|
||||
|
||||
Roberts, M., Hayes, W., Hunt, B.R., Mount, S.M. & Yorke, J.A. (2004). <a href="https://doi.org/10.1093/bioinformatics/bth408">Reducing storage requirements for biological sequence comparison</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 20, 3363–3369.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Stafford2011" class="csl-entry">
|
||||
|
||||
Stafford, D. (2011). <span style="font-style: italic;"><a href="https://zimbry.blogspot.com/2011/09/better-bit-mixing-improving-on.html">Better Bit Mixing - Improving on MurmurHash3's 64-bit Finalizer</a></span>. Available at: <a href="https://zimbry.blogspot.com/2011/09/better-bit-mixing-improving-on.html"><https://zimbry.blogspot.com/2011/09/better-bit-mixing-improving-on.html></a>. Last accessed 12 September 2026.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Steele2014-cz" class="csl-entry">
|
||||
|
||||
Steele, G.L., Jr, Lea, D. & Flood, C.H. (2014). <a href="https://doi.org/10.1145/2660193.2660195">Fast splittable pseudorandom number generators</a>. In: <span style="font-style: italic;">Proceedings of the 2014 ACM International Conference on Object Oriented Programming Systems Languages & Applications</span>. ACM, New York, NY, USA.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Zheng2020-ji" class="csl-entry">
|
||||
|
||||
Zheng, H., Kingsford, C. & Marçais, G. (2020). <a href="https://doi.org/10.1093/bioinformatics/btaa472">Improved design and analysis of practical minimizers</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 36, i119–i127.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Zheng2021-cc" class="csl-entry">
|
||||
|
||||
Zheng, H., Kingsford, C. & Marçais, G. (2021). <a href="https://doi.org/10.1093/bioinformatics/btab313">Sequence-specific minimizers via polar sets</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 37, i187–i195.
|
||||
|
||||
</div>
|
||||
|
||||
</div>
|
||||
@@ -0,0 +1,33 @@
|
||||
# Super-kmers
|
||||
|
||||
A **super-kmer** is a maximal run of consecutive, overlapping kmers from a read that share the same canonical minimizer (see [Minimizer selection](theory-kmer_indexing-minimizer_selection)). Each kmer in the run overlaps the next by $k-1$ nucleotides. A super-kmer is capped at 256 nucleotides; a longer run is split at that boundary.
|
||||
|
||||
For a random minimizer of length $m$ over kmers of length $k$, the expected length of a super-kmer is approximately ([Golan & Shur 2025](#ref-Golan2025-xf); [Zheng et al. 2020](#ref-Zheng2020-ji)):
|
||||
|
||||
$$L_{\text{nt}} \approx \frac{k-m+2}{2} + k - 1$$
|
||||
|
||||
For $k=31$, $m=13$ this is about 40 nucleotides; in practice super-kmers rarely exceed a few dozen nucleotides.
|
||||
|
||||
## Canonical super-kmers
|
||||
|
||||
A **canonical super-kmer** is the canonical form of a super-kmer (see [DNA encoding](theory-kmer_indexing-encoding)): the lexicographic minimum of the super-kmer and its reverse complement. When a read and its reverse complement are both encountered, they produce super-kmers that are reverse complements of each other; both reduce to the same canonical super-kmer, so a genomic region is represented once regardless of which strand was read.
|
||||
|
||||
Super-kmers are the unit of work used throughout construction and querying: sequences are decomposed into super-kmers first, and every downstream step (partition routing, deduplication, counting) operates on them rather than on individual kmers.
|
||||
|
||||
## Bibliography
|
||||
|
||||
<div id="refs" class="references csl-bib-body hanging-indent">
|
||||
|
||||
<div id="ref-Golan2025-xf" class="csl-entry">
|
||||
|
||||
Golan, S. & Shur, A.M. (2025). <a href="https://doi.org/10.1007/978-3-031-82670-2\_25">Expected density of random minimizers</a>. In: <span style="font-style: italic;">Lecture Notes in Computer Science</span>, Lecture Notes in Computer Science. Springer Nature Switzerland, Cham, pp. 347–360.
|
||||
|
||||
</div>
|
||||
|
||||
<div id="ref-Zheng2020-ji" class="csl-entry">
|
||||
|
||||
Zheng, H., Kingsford, C. & Marçais, G. (2020). <a href="https://doi.org/10.1093/bioinformatics/btaa472">Improved design and analysis of practical minimizers</a>. <span style="font-style: italic;">Bioinformatics (Oxford, England)</span>, 36, i119–i127.
|
||||
|
||||
</div>
|
||||
|
||||
</div>
|
||||
@@ -0,0 +1,19 @@
|
||||
# Whole-index distance metrics
|
||||
|
||||
These metrics operate on the kmer sets or counts stored in an index directly — no per-locus SNP calling, no sibling annex required. `A`/`B` denote the two genomes being compared; $c_i^A$/$c_i^B$ are their raw counts at kmer $i$, $p_i^A$/$p_i^B$ the corresponding relative frequencies ($p_i = c_i / \sum_j c_j$).
|
||||
|
||||
| Metric | 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 relative frequencies $p_i$ 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.
|
||||
|
||||
See [`phylo`](usage-phylo) for how to select a metric and the output formats.
|
||||
@@ -0,0 +1,47 @@
|
||||
# Central-position SNP model
|
||||
|
||||
## 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.
|
||||
|
||||
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.
|
||||
|
||||
## Substitution corrections (`snp-*`)
|
||||
|
||||
For a genome pair, let $L$ be its total number of eligible loci, $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$$
|
||||
|
||||
A rate-heterogeneity correction applies to every formula above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (same weight, same $x$) — see [Gamma rate heterogeneity](theory-phylogeny-snp_based-gamma_rate_heterogeneity) for the correction itself and how $\alpha$ is estimated.
|
||||
|
||||
See [`phylo`](usage-phylo) for how to select a `snp-*` value and the prerequisite index annex it requires.
|
||||
@@ -0,0 +1,27 @@
|
||||
# Gamma rate heterogeneity
|
||||
|
||||
Real substitution rates vary across sites rather than being uniform, which a plain substitution correction (see [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model)) ignores. The standard $+\Gamma$ correction models this by replacing each $-\ln(x)$ term in a correction formula with $\alpha\left(x^{-1/\alpha}-1\right)$ — the same $x$, the same weight, under the assumption that the per-site rate is $\mathrm{Gamma}(\alpha,\alpha)$-distributed (mean 1) rather than fixed. A small $\alpha$ indicates strong among-site rate heterogeneity (many near-invariant sites, a few fast ones); as $\alpha$ grows large the correction converges to the uncorrected formula.
|
||||
|
||||
## Automatic $\alpha$ estimation
|
||||
|
||||
$\alpha$ can be estimated automatically from the index itself, 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 in the correction itself, 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.
|
||||
|
||||
**This is not Jin & Nei's (1990) original estimator.** Their publication defines the $+\Gamma$ distance formula itself and, absent an estimate, recommends the fixed default $\alpha = 1$ — it does not propose a way to estimate $\alpha$ from data. The method-of-moments estimator above is a standard consequence of the Poisson–Gamma relationship between substitution counts and gamma-distributed rate variation, applied per-partition; it is `obikmer`'s own addition, not part of the cited correction.
|
||||
|
||||
See [`phylo`](usage-phylo#distance-matrix---distance) for how to select a fixed $\alpha$ or request automatic estimation.
|
||||
@@ -0,0 +1,24 @@
|
||||
# Sampling design
|
||||
|
||||
Exhaustive computation over every non-monomorphic family in an index (see [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model)) is the default, but can be bounded to approximately `N` families instead. Sampling is designed around two properties: staying representative of the whole index, and letting the user favor families that actually carry signal.
|
||||
|
||||
## Proportional sampling
|
||||
|
||||
Families are drawn in proportion to how many candidate families each part of the index actually holds, rather than uniformly across index partitions — this keeps the sample's composition representative of the whole index regardless of how unevenly candidate families happen to be distributed across partitions.
|
||||
|
||||
## Per-family informativeness: Shannon entropy
|
||||
|
||||
A family's informativeness is measured directly rather than assumed uniform. For a family, entropy is computed over its observed states across genomes:
|
||||
|
||||
| Quantity | Definition |
|
||||
| --- | --- |
|
||||
| `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` | the same, reduced to the 4 plain bases |
|
||||
|
||||
A family where every genome carries the same single state has entropy 0 (uninformative); one where genomes are spread evenly across several states has higher entropy (more informative for distinguishing genomes).
|
||||
|
||||
## Entropy-biased sampling
|
||||
|
||||
By default, sampling draws families uniformly. Weighting instead by informativeness biases the draw toward a target entropy $\mu$ with a Gaussian curve of width $\sigma$: a family's probability of being kept scales with $\exp\!\left(-\dfrac{(\text{entropy15} - \mu)^2}{2\sigma^2}\right)$ — a soft preference, not a hard cutoff, so no family is categorically excluded purely for having low or high entropy.
|
||||
|
||||
See [`phylo`](usage-phylo#sampling-at-scale---subsample---shannon---entropy) for how to bound and bias the sample, and for session/checkpointing behavior.
|
||||
@@ -0,0 +1,22 @@
|
||||
# The 16-state Sankoff model
|
||||
|
||||
Each family (see [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model)) is treated as a character with 16 possible states: one per subset of the 4 possible central bases, including the empty subset (no observed form). A genome pair's calibration reduces to two separately-estimated first-order Markov models, both restricted to genome pairs below a configurable SNP-ratio ceiling (excluding saturated pairs from skewing the calibration):
|
||||
|
||||
- **Cardinality transitions** ($5 \times 5$, row-stochastic): the probability that a family observed with $i$ forms in one genome ($i \in \{0,\dots,4\}$) is observed with $j$ forms in the other.
|
||||
- **Composition transitions** ($4 \times 4$, row-stochastic): the probability that an unambiguous, single-copy base observed as one of A/C/G/T in one genome is observed as another base in the other genome.
|
||||
|
||||
## Combining into a 16×16 cost matrix
|
||||
|
||||
The two matrices are *not* combined as a plain tensor product. For a pair of states $(A, B)$ — each a subset of $\{A,C,G,T\}$ — let `shared` $= A \cap B$, `lost` $= A \setminus B$, `gained` $= B \setminus A$. The log-probability of the transition $A \to B$ is:
|
||||
|
||||
$$\log P(A,B) = \underbrace{\log P_{\text{card}}(|A|,|B|)}_{\text{cardinality term, omitted when loss events are treated as free}} \;+\; \underbrace{\sum_{i \,\in\, \text{shared}} \log P_{\text{comp}}(i,i)}_{\text{bases conserved on both sides}} \;-\; \underbrace{\text{best\_pairing\_cost}(\text{lost}, \text{gained})}_{\text{substitutions among the differing bases}}$$
|
||||
|
||||
`best_pairing_cost` enumerates every injection pairing elements of `lost` with elements of `gained` and keeps the one minimizing $-\sum \log P_{\text{comp}}$ over the pairs — the discrete analogue of "prefer a substitution over an independent loss and gain": if a base is lost from one side and a different base is gained on the other, that is scored as a single substitution between them (cheaper, under any sane calibration, than treating them as two unrelated cardinality-changing events) whenever such a pairing is possible.
|
||||
|
||||
The resulting $16 \times 16$ log-probability matrix is row-normalized in log-space (log-sum-exp, not a direct `exp()`/divide, for numerical stability), giving a genuine row-stochastic transition matrix $P$. The cost matrix is then $\text{cost}(A,B) = -\ln P(A,B)$, symmetrized by simple averaging:
|
||||
|
||||
$$\text{sym}(A,B) = \frac{\text{cost}(A,B) + \text{cost}(B,A)}{2}$$
|
||||
|
||||
This calibrated cost matrix is what every downstream export (for TNT, PhyG, IQ-TREE) reuses, recoded into each tool's own format.
|
||||
|
||||
See [`phylo`](usage-phylo#sankoff-calibration-and-phylogenetic-exports) for how to run the calibration and produce these exports.
|
||||
+7
-85
@@ -17,7 +17,7 @@ obikmer phylo INDEX [OPTIONS]
|
||||
| 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 |
|
||||
| `--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 |
|
||||
@@ -25,85 +25,9 @@ obikmer phylo INDEX [OPTIONS]
|
||||
| `--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:
|
||||
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).
|
||||
|
||||
### 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.
|
||||
`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).
|
||||
|
||||
### Output
|
||||
|
||||
@@ -120,7 +44,7 @@ Neighbor-Joining and UPGMA trees (`--nj`/`--upgma`) are always built from every
|
||||
|
||||
## Central-position SNP model
|
||||
|
||||
Requires the sibling annex, built once per index:
|
||||
See [Central-position SNP model](theory-phylogeny-snp_based-central_snp_model) for the family/eligibility definitions. Requires the sibling annex, built once per index:
|
||||
|
||||
| Option | Description |
|
||||
| --- | --- |
|
||||
@@ -131,8 +55,6 @@ Requires the sibling annex, built once per index:
|
||||
| `--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).
|
||||
@@ -181,7 +103,7 @@ A `--sankoff`-family run and a `snp-*` `--distance` run share the same cached sa
|
||||
|
||||
If a run using `--session` is interrupted (crash, kill, `Ctrl-C`), the next run against the same `--session DIR` resumes from the last automatic checkpoint (roughly every 8 partitions'-worth of sampling progress) instead of starting over. The resumed run's sample is **not** guaranteed to be identical to what an uninterrupted run would have produced past that checkpoint — each process draws its own independent random sequence, same as any two separate invocations do — but nothing already checkpointed is lost, and no work needs redoing beyond that point.
|
||||
|
||||
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.
|
||||
Without `--subsample`, every variable family (family size ≥ 2) is used. With `--subsample N`, roughly `N` families are kept instead (see [Sampling design](theory-phylogeny-snp_based-sampling_design) for how the draw 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
|
||||
|
||||
@@ -200,7 +122,7 @@ Run with `--subsample N --shannon` to get a bounded diagnostic sample instead of
|
||||
|
||||
### `--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).
|
||||
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 (see [Sampling design](theory-phylogeny-snp_based-sampling_design) for the weighting formula). 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.
|
||||
|
||||
@@ -219,7 +141,7 @@ The first `phylo` run on a given index that uses `--entropy`/`--entropy-sd` pays
|
||||
|
||||
### 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)$.
|
||||
Each family is a character with 16 possible states: one per subset of the 4 possible central bases (including the empty subset). See [The 16-state Sankoff model](theory-phylogeny-snp_based-sankoff_model) for how the calibrated cost matrix is built from cardinality and base-composition transitions.
|
||||
|
||||
`--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.
|
||||
|
||||
|
||||
Reference in new issue
Block a user