obikmer/docmd/theory/entropy.md

Kmer entropy filter

Low-complexity kmers (polyA, polyT, tandem repeats) are detected and excluded during phase 1. The filter computes a normalized Shannon entropy over sub-words of multiple sizes, corrected for one source of bias: the small number of observations within a single kmer relative to the number of possible sub-words.

Sub-word frequencies

For a kmer of length k and a sub-word size ws (1 ≤ ws ≤ ws_max, typically ws_max = 6), extract the n_{\text{words}} = k - ws + 1 overlapping sub-words by sliding a window of length ws:

w_i = \text{kmer}[i \mathinner{..} i+ws-1], \quad i = 0, \ldots, n_{\text{words}}-1

Each sub-word is tallied under its own raw 2-bit-packed value — no canonicalization. Let f_j be the count of raw word j among the n_{\text{words}} sub-words (\sum_j f_j = n_{\text{words}}), over the 4^{ws} possible raw words.

An earlier version of this filter first folded each sub-word into a circular+reverse-complement equivalence class, then "unfolded" the observed class frequency back onto its members to correct for unequal class sizes. That machinery bought nothing it was claimed for — see Why no equivalence classes below — while measurably weakening detection of the very sequences the filter exists to catch, so it was removed.

Corrected Shannon entropy

H_{\text{corr}} = \log(n_{\text{words}}) - \frac{1}{n_{\text{words}}} \sum_j f_j \log f_j

This is a plain Shannon entropy over the observed raw-word frequencies.

Maximum entropy correction for small samples

With only n_{\text{words}} observations over 4^{ws} possible raw words, the achievable maximum entropy is bounded by the most uniform integer distribution over 4^{ws} categories.

Let c = \lfloor n_{\text{words}} / 4^{ws} \rfloor and r = n_{\text{words}} \bmod 4^{ws}. The most uniform integer distribution assigns frequency c+1 to r categories and c to the remaining 4^{ws} - r, with the convention 0 \log 0 = 0:

H_{\max} = -\left[(4^{ws} - r)\,\frac{c}{n_{\text{words}}}\log\frac{c}{n_{\text{words}}} + r\,\frac{c+1}{n_{\text{words}}}\log\frac{c+1}{n_{\text{words}}}\right]

When n_{\text{words}} < 4^{ws}: c=0, r=n_{\text{words}}, and the formula reduces to H_{\max} = \log(n_{\text{words}}) — a single unified expression covers both regimes. A truly random sequence achieves H_{\text{corr}} \approx H_{\max}.

Normalized entropy

\hat{H}(ws) = \frac{H_{\text{corr}}}{H_{\max}} \in [0, 1]

Final score

The filter computes \hat{H}(ws) for each word size ws from 1 to ws_max and returns the minimum:

\text{entropy}(kmer) = \min_{ws=1}^{ws_{\max}} \hat{H}(ws)

A value near 0 indicates low complexity (e.g. AAAA…); near 1 indicates high complexity. A kmer is rejected if \text{entropy}(kmer) < \theta, where \theta is a collection parameter (default 0.7). The minimum across word sizes ensures that any scale of repetition is detected independently: polyA is caught at ws=1, dinucleotide repeats at ws=2, etc.

Why no equivalence classes

A prior design folded each sub-word into the canonical form of its circular-rotation + reverse-complement equivalence class before tallying, on the reasoning that (a) it guarantees \text{entropy}(K) = \text{entropy}(\text{revcomp}(K)), and (b) collapsing phase-shifted repeats (e.g. ATG ≡ TGA ≡ GAT) into one class better reflects that they are "the same" low-complexity pattern.

Both properties already hold for the raw, unfolded entropy above, without any class machinery:

• Reverse complement: for any K of length n, window j of \text{revcomp}(K) equals \text{revcomp} of window (n{-}ws{-}j) of K. This is a bijection between the window sets under which each window maps to its own revcomp — and revcomp is itself a bijection (involution) on the space of raw ws-mers. So the multiset of raw-word frequencies for \text{revcomp}(K) is exactly a relabeling of the multiset for K, and Shannon entropy — a function of the frequency multiset alone — is exactly invariant. No folding required, for any K.
• Tandem repeats: a period-p repeat sampled by a stride-1 sliding window naturally cycles through its own rotations as raw tokens (e.g. ATGATGATG… yields the raw words ATG, TGA, GAT in rotation as the window slides). The low diversity this represents (few distinct raw words out of 4^{ws} possible) is already visible in the raw frequency distribution — no folding needed to detect it.

What the fold-then-unfold step actually did was credit each observed class with the frequency of equivalence-class members that were never observed on the read strand, inflating H_{\text{corr}} for genuine repeats. Worked example: k=31, ws=3, kmer = ATG repeated (n_{\text{words}}=29, all 29 windows fall into one class of size 6 under the old scheme — 3 rotations × forward/revcomp):

	H_{\text{corr}}	normalized
old (folded, class size 6)	\log 6 \approx 1.79	\approx 0.53
current (raw, unfolded)	\log 3 \approx 1.10	\approx 0.33

The gap is not a rounding artifact: per sub-word order, the folded score for this same repeat swings from 0.53 (ws=3, aligned with the period) up to 1.03 (ws=5, misaligned with the period) — i.e. a period-3 repeat could score above the theoretical maximum for a random sequence, depending on which ws happens to divide the repeat's period. The raw formula stays flat at ≈0.33–0.40 across ws=2..6 regardless of alignment, which is the robustness the "minimum across ws" design was meant to provide in the first place.

Interpretation as an effective number of classes

H_{\text{corr}} is a standard Shannon entropy over raw words, so the classical perplexity interpretation holds directly: N_{\text{eff}} = e^{H_{\text{corr}}} is the number of equiprobable raw words that would yield the same entropy.

For the normalised score \hat{H}, dividing by H_{\max} changes the logarithm base:

\hat{H} = \frac{\log N_{\text{eff}}}{\log N_{\max}} = \log_{N_{\max}} N_{\text{eff}} \quad \Longleftrightarrow \quad N_{\text{eff}} = N_{\max}^{\,\hat{H}}

The property is preserved: \hat{H} is the logarithm (in base N_{\max}) of the effective number of equi-represented raw words.

In the large-sample limit (n_{\text{words}} \gg 4^{ws}), N_{\max} \approx 4^{ws}, giving:

N_{\text{eff}} \approx 4^{ws \cdot \hat{H}}

This has a clean interpretation: ws \cdot \hat{H} is the effective word length (in bases) of a perfectly uniform distribution that would produce the same entropy. At \hat{H} = 1 the full space of 4^{ws} words is used; at \hat{H} = 0.5 with ws=2, only 4^1 = 4 effective words out of 16 are occupied.

In our actual regime, n_{\text{words}} is small and 4^{ws} can exceed n_{\text{words}}, so H_{\max} < \log(4^{ws}) due to the small-sample correction. The exact effective count is N_{\max}^{\hat{H}}, not 4^{ws \cdot \hat{H}}.

Properties

The entropy score is a function of the kmer sequence alone — it does not depend on the surrounding context or on the position within any genome. Two consequences:

• Orientation invariance: \text{entropy}(K) = \text{entropy}(\text{revcomp}(K)) — see Why no equivalence classes above for why this holds without any explicit strand-folding step.
• Context independence: the same kmer is always rejected or always kept, regardless of which genome it occurs in, where in that genome it appears, or which strand is considered. The filter defines a fixed partition of the kmer space into low-complexity and valid kmers.