Characterising the D2 statistic: word matches in biological sequences

Word matches are often used in sequence comparison methods, either as a measure of sequence similarity or in the first search steps of algorithms such as BLAST or BLAT. The D2 statistic is the number of matches of words of k letters between two sequences. Recent advances have been made in the characterisation of this statistic and in the approximation of its distribution. Here, these results are extended to the case of approximate word matches. We compute the exact value of the variance of the D2 statistic for the case of a uniform letter distribution, and introduce a method to provide accurate approximations of the variance in the remaining cases. This enables the distribution of D2 to be approximated for typical situations arising in biological research. We apply these results to the identification of cis-regulatory modules, and show that this method detects such sequences with a high accuracy. The ability to approximate the distribution of D2 for both exact and approximate word matches will enable the use of this statistic in a more precise manner for sequence comparison, database searches, and identification of transcription factor binding sites.

💡 Research Summary

The paper provides a comprehensive theoretical and practical treatment of the D₂ statistic, a widely used alignment‑free measure of similarity between biological sequences that counts the number of matching k‑mers (or approximate matches allowing up to t mismatches) shared by two sequences. While earlier work derived the mean of D₂ and offered asymptotic approximations for its variance and distribution (often assuming a normal or gamma law), those results were limited to exact matches, uniform nucleotide composition, and did not address the distribution’s tail behavior, which is crucial for significance testing.

The authors first formalize D₂(n_A, n_B, k, t, η) where n_A and n_B are the lengths of the two sequences, k is the word length, t is the allowed number of mismatches, and η parameterizes a strand‑symmetric Bernoulli model (ξ_A = ξ_T = (1+η)/4, ξ_G = ξ_C = (1−η)/4). They impose periodic boundary conditions to simplify calculations, noting that the same results apply to linear sequences after appropriate preprocessing.

Mean. By extending the perturb‑binomial framework, they derive a closed‑form expression for the expected value of D₂ that holds for any t (0 ≤ t ≤ k) and any η. The formula reduces to the classic n_A n_B 4⁻ᵏ C(k, t) term when η = 0 and t = 0, but more generally incorporates the bias in nucleotide frequencies.

Variance. The variance is expressed as a double sum over all pairs of word positions, which they decompose into four dependency neighborhoods: the “accordion” set (overlap in both sequences), further split into diagonal and off‑diagonal components, and the “crabgrass” set (overlap in only one sequence). For uniform letters (η = 0) and exact matches (t = 0) the variance collapses to known results (e.g., Eq. 20 of Reinert et al., 2009). For non‑uniform letters and exact matches they recover the expression from Waterman et al. (2010).

When t > 0 (approximate matches) the off‑diagonal accordion term becomes analytically intractable for η ≠ 0. The authors therefore introduce a function Φ(k, t, η) that captures the entire accordion contribution. Φ is estimated by massive Monte‑Carlo simulations of paired sequences of length 2k − 1 for k ≤ 16 and all admissible t, generating 10⁹ pairs to obtain high‑precision tables. They show that for n_A = n_B = 2k − 1 the variance simplifies to n_A n_B Φ(k, t, η), enabling rapid variance approximation for realistic sequence lengths.

Distribution Approximation. Prior work suggested a normal approximation for large n, but the authors demonstrate that for biologically relevant sizes (n = 100–3200, k = 2–16) the normal or even gamma approximations mis‑represent the extreme right tail, leading to inaccurate p‑values. They propose fitting a scaled beta distribution B(α, β) on the interval