The distribution of word matches between Markovian sequences with periodic boundary conditions

Abstract

Word match counts have traditionally been proposed as an alignment-free measure of similarity for biological sequences. The D(2) statistic, which simply counts the number of exact word matches between two sequences, is a useful test bed for developing rigorous mathematical results, which can then be extended to more biologically useful measures. The distributional properties of the D(2) statistic under the null hypothesis of identically and independently distributed letters have been studied extensively, but no comprehensive study of the D(2) distribution for biologically more realistic higher-order Markovian sequences exists. Here we derive exact formulas for the mean and variance of the D(2) statistic for Markovian sequences of any order, and demonstrate through Monte Carlo simulations that the entire distribution is accurately characterized by a Pólya-Aeppli distribution for sequence lengths of biological interest. The approach is novel in that Markovian dependency is defined for sequences with periodic boundary conditions, and this enables exact analytic formulas for the mean and variance to be derived. We also carry out a preliminary comparison between the approximate D(2) distribution computed with the theoretical mean and variance under a Markovian hypothesis and an empirical D(2) distribution from the human genome.

FIG. 1.

Covering of the sequence X with θ-mers for the calculation of (a) in the case where k ≥ θ, and (b) in the case where k < θ.

FIG. 2.

Contributions to Var(D₂) via the sum in Equation (26). The left-hand diagram shows the (i′, j′)-plane for a fixed value of (i, j), shown as the black square. The right-hand diagram is an expanded view of the “accordion” region −k + 1 ≤ s, t ≤ k − 1, where t = i′ − i and s = j′ − j up to PBCs [see Equations (41) and (42)].

FIG. 3.

The exact distribution of the D₂ statistic for short sequences of length m, n and words of length k from a Markov model of order θ and alphabet of size d. The Markov matrix M has been generated randomly in each case, and the exact distribution has been calculated by enumerating all d^m+n possible sequence pairs. Also shown (dashed curve) is the cumulative distribution of the Pólya-Aeppli distribution with mean and variance set to the theoretical values using the formulas of Section 3.

FIG. 4.

Comparison of empirical cumulative distribution function for simulated D₂ using Markov matrices for human chromosome 1 from Chor et al. (2009a,b), for orders θ = 0, 1, 2, and 3. 10,000 pairs per order, word length k = 8, alphabet size d = 4, sequence lengths m = n =1,000.

FIG. 5.

Comparison of Pólya-Aeppli, gamma, and normal cumulative distributions with empirical cumulative distribution function for simulated D₂ using Markov order 3 matrix for human chromosome 1 from Chor et al. (2009a,b). 10,000 pairs, word length k = 8, sequence lengths m = n = 1,000.

FIG. 6.

Density of the empirical distribution of D₂ from human chromosome 1 sample data from Ensembl compared with gamma distributions with calculated mean and variance, based on Markov models of various orders θ. 10,000 sample pairs, word length k = 8, and sequences lengths m = n = 1,000 (upper plot), word length k = 5 and sequences lengths m = n = 300 (lower plot).

FIG. 7.

Arrangements of word matches contributing to (a) nonoverlapping words, V₀, (b) crabgrass V₁, (c) accordions V₂, V₃, and (d) accordion V₄. Images of words due to periodic boundary conditions are shown as a dashed outline.

FIG. 8.

Arrangements of word matches contributing to subcase 3(i) when ρ = (k − s) mod (t − s)> 0 (upper figure) and ρ = 0 (lower figure).

FIG. 9.

Arrangements of word matches contributing to V4 (first four cases).

FIG. 10.

Arrangements of word matches contributing to V4 (final three cases).