Table of Contents
Gamma rate heterogeneity
Real substitution rates vary across sites rather than being uniform, which a plain substitution correction (see Central-position 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 for how to select a fixed \alpha or request automatic estimation.
Wiki sidebar
Theory
Kmer indexing
- DNA encoding
- Kmers
- Minimizer selection
- Super-kmers
- Partitioning and indexing architecture
- Low-complexity kmer filter
Phylogeny
Kmer-based
SNP-based
Usage
- superkmer
- index
- merge
- filter
- select
- query
- dump
- annotate
- phylo
- unitig
- estimate
- convert
- utils
- pack
- Predicates and taxonomy paths