Estimation of Substitution and Indel Rates via $k$-mer Statistics

Estimating the mutation rate between two evolutionarily related sequences is a classical question in molecular evolution, with roots that pre-date the genomics era [20]. Early quantitative efforts focused on amino-acid substitution: the seminal PAM matrices of Dayhoff et al. converted curated alignments of close homologues into an evolutionary time-scale [3], while the BLOSUM series by Henikoff and Henikoff mined ungapped blocks of conserved proteins to improve sensitivity for more diverged sequences [6]. These approaches and later profile-based HMM models [4] were derived from pairwise or multiple alignments and remain the gold standard when accurate alignments are available.

Over the last decade, however, high-throughput sequencing has shifted the scale of comparative genomics from dozens to millions of genomes, rendering high computational cost pipelines (e.g. quadratic-time) increasingly impractical. Consequently, alignment-free techniques that summarize sequences by inexpensive statistics have become indispensable [17, 21]. These approaches most commonly utilize $k$-mer sets and sketches thereof. Popular tools such as Mash [10], Skmer [13] and more recent sketch-corrected frameworks like Sylph [16] and FracMinHash-based methods [7, 8, 11, 15] can build whole-genome phylogenies, screen metagenomic samples, and compute millions of pairwise point-mutation rate estimates in minutes rather than days.

Despite their empirical success, theoretical understanding of alignment-free estimators has lagged behind practice. Nearly all existing models treat evolution as a pure-substitution process, ignoring insertions and deletions (indels), or else their performance in the presence of indels is often not thoroughly evaluated [12]. When indels are frequent, substitution-only estimators systematically inflate divergence and can misplace taxa – even in otherwise well-resolved trees of primates constructed from $k$-mer Jaccard similarities [10]. Recent work has quantified how $k$-mer-based statistics are affected by substitutions and are also used to estimate substitution-only mutation rates [8], yet a principled treatment that jointly infers substitution and indel parameters from $k$-mer statistics is still absent. This omission is particularly significant because indels represent a substantial fraction of genomic variation and play crucial roles in evolution [19]. Such indels cause substitution-only $k$-mer-based methods to underperform, as just like with substitutions, disruption of $k$-mer content by indel events affects at least $k$ $k$-mers, often leading to overestimates of mutation rates [2, 8].

In this paper we introduce the first closed-form, alignment-free estimators for the three fundamental mutation parameters: substitution rate $p_s$, deletion rate $p_d$, and mean insertion length $d$ under a model that explicitly incorporates single-nucleotide substitutions, deletions, and geometrically-distributed insertions. Starting from elementary counts of unmutated and single-deletion $k$-mers, we derive algebraic expressions for $p_s$, $p_d$, and $d$ and prove a sub-exponential concentration bound that guarantees a strong form of consistency as the sequence length grows as detailed in our main contribution: Theorem 9. Simulations on synthetic and real bacterial genomes demonstrate that modeling indels yields markedly more accurate distance estimates than substitution-only approaches. The remainder of this paper is organized as follows: Section 2 defines our mutation model, Section 3 derives basic statistical properties, Section 4 presents our estimators, Section 5 provides theoretical guarantees about these estimators, Section 6 details our implementation, Section 7 presents experimental results, and Section 8 concludes with a discussion of implications and future directions.

We first describe the mutation model under consideration. We use the model from [14]. This model is a type of indel channel and various variations of it have been used to model sequence evolution (e.g. [5]). Let $S$ be a string over the alphabet $\{A,C,G,T\}$, and let $L$ be the number of characters in $S$. Let $S_i$ denote the $i$-th character in $S$ where $1\le i\le L$. We define four operations on a string $S$ as follows:

• $\mathrm{■}$

  Substitute($i$): select a character $c\in \{A,C,G,T\}\setminus \{S_i\}$ uniformly at random and replace $S_i$ by $c$.

• $\mathrm{■}$

  Delete($i$): remove the character $S_i$.

• $\mathrm{■}$

  Stay($i$): do nothing.

• $\mathrm{■}$

  Insert($i$, $s$): insert string $s$ between $S_{i-1}$ and $S_i$ if $i>0$. If $i=0$, prepend $s$ to the start of $S$ instead.

The mutation process takes as input a string $S$ and three parameters $p_s$, $p_d$, and $d$, where , , and $d\ge 0$. Then,

1. 1.

   For each $i$, let $a_i$ be the operation to be performed at position $i$. Then, $a_i=$Sub with probability $p_s$, $a_i=$Del with probability $p_d$, or $a_i=$Stay with the remaining probability $1-p_s-p_d$. (These operations are not performed at this point, but only recorded.)

2. 2.

   Let $track$ be a function mapping from a position in the original string $S$ to its position in the modified $S$. Initially, $track(i)=i$ for all $i$. We assume that $track$ is updated accordingly whenever an insertion or deletion operation is performed.

3. 3.

   For each $i$ that is a substitute action, apply Substitute($i$).

4. 4.

   For each $i$, let $l_i\ge 0$ be a sample from a geometric distribution with mean $d$. Then generate a random string $Q_i$ of length exactly $l_i$ by drawing each character in $Q_i$ from $\{A,C,G,T\}$ independently and uniformly at random. This is equivalent to sampling a uniformly random string among all strings of length $l_i$ with characters in $\{A,C,G,T\}$. If $l_i=0$, $Q_i$ is the empty string.

5. 5.

   For every $i$, execute Insert($track(i)$, $Q_i$).

6. 6.

   For every $i$ such that $a_i=$ Del, execute Delete($track(i)$).

7. 7.

   Return the resulting $S^{\prime }$.

In this section, we define several quantities whose statistics we use to design an estimator of the parameters of the mutation process. For the remainder of the manuscript, we assume that the mutation process described in Section 2 is applied to string $S$ of length $L$, with the unknown parameters $p_s$, $p_d$, and $d$, and a mutated string $S^{\prime }$ is returned. The theoretical results we present are centered around the concept of $k$-spans. We define ${𝒦}_i$, the $k$-span at position $1\le i\le L-k+1$ as the range of integers $[i,i+k-1]$ (inclusive of the endpoints of the range). In simpler terms, a $k$-span captures the interval of a $k$-mer. We assume that the string $S$ has at least $k$ nucleotides, and therefore, at least one $k$-mer (and at least one $k$-span). For the sake of theoretical rigor, we use $k$-spans to develop our results, and later discuss how our practical implementation uses $k$-mers to estimate the substitution and indel rates.

We can analogously obtain expressions for the expectations of the corresponding counts of “C”s, “G”s, and “T”s in $S^{\prime }$, denoted by $f_C^{}{}_{}{}^{\prime }$, $f_G^{}{}_{}{}^{\prime }$, and $f_T^{}{}_{}{}^{\prime }$, respectively.

Next, we say that a $k$-span has no mutations if (a) for all $k$ positions in the $k$-span, the mutation process picks a Stay($\cdot$) action, and if (b) for all $k-1$ intermediate positions in the $k$-span, the mutation process does not insert anything. Let $𝒩$ be the number of such $k$-spans. Let $𝒦$ be the number of all $k$-spans in $S$.

The last quantities we will use correspond to the number of $k$-spans with a single kind of mutation. Let $𝒮$ be the number of $k$-spans with a single substitution, and no other mutations, let $𝒟$ be the number of $k$-spans with a single deletion and no other mutations, and let $ℐ$ be the number of $k$-spans with a single insertion and no other mutations.

We describe next our estimators for the parameters of the mutation model $p_s$, $p_d$, and $d$. We derive estimators for these rates based on the statistics defined in Section 3. Recall that $S$ is the known input string to the mutation process, and $S^{\prime }$ is the resulting random string or observation.

Our mutation model assumes , and $d\ge 0$. In addition, we assume that $L\ge k$ which implies that $𝐄[𝒩]>0$. In reality, strings are typically much longer than $k$-mer sizes and therefore this is a reasonable assumption.

From Lemmas 3 and 4 we get

From (1), (2) and (6), we obtain the following system of linear equations with variables $p_s$, $p_d$, and $d$.

Solving this system of equations, we obtain that:

Given $S$ and $S^{\prime }$, we can compute $L^{\prime }$ and $f_A^{}{}_{}{}^{\prime }$. Let us also assume that we know $𝒩$ and $𝒟$ (we will discuss how to find these in Section 6). By replacing the expectations above with these observations, we obtain our estimators for the parameters of the mutation model. That is,

We briefly comment on our choice of estimators for the mutation model parameters, as various statistical approaches based on a different set of observables could yield a different set of estimators. For example, we considered a variant of the estimators based on the counts of $k$-spans with a single insertion, deletion, or substitution (i.e., $ℐ$, $𝒮$, and $𝒟$). These quantities contain enough information to estimate the mutation parameters, specifically, by solving the non-linear system of equations given by (3), (4), and (5); the resulting estimators performed quite well in real data. However, establishing theoretical guarantees for these estimators proved challenging, as they were defined as roots of degree-$k$ polynomials . Our current estimators address this theoretical limitation as they involve solely linear equations. As we shall see in Section 7, the performance of our estimators in real data is strong, and they strike a more favorable balance by offering both reasonable accuracy and rigorous theoretical guarantees.

In this section, we provide a theoretical guarantee for the estimator (10) of the substitution rate. In particular, we show that our estimate of $p_s$ is not only a consistent estimator, but also that it is tightly concentrated in a symmetric interval around the true value of $p_s$. Similar techniques can likely be used to prove the consistency and concentration of $p_d$ and $d$, though we do not do so in this paper. We provide the bias analysis for $\widehat{p_s}$ by proving asymptotically tight concentration bounds for $𝒩$, $𝒟$, and $L^{\prime }-4f_A^{}{}_{}{}^{\prime }$.

The concentration bounds for $𝒩$ and $𝒟$ stem from the fact that, from the perspective of the mutation process $k$-spans located more than $k$ positions apart are independent of each other.

We prove a similar result for $P=L^{\prime }-4f_A^{}{}_{}{}^{\prime }$ by noting that $P$ can expressed as a sum of independent sub-exponential random variables. Our concentration bound for $P$ is as follows.

Finally, we prove that $L^{\prime }$ is also strongly concentrated around $𝐄[L^{\prime }]$.

We can piece the results of these lemmas together and prove the following result.

The requirement that in this theorem does not restrict generality: aside from equal nucleotide frequency (where the estimators are naturally invalid), at least one character $c\in \{A,C,T,G\}$ must satisfy $4f_c\le L$. In addition, the assumption that holds when $p_s$ and $p_d$ are small (e.g., and suffices) which is the case most frequently encountered in practice.

Theorem 9 is our central theoretical result, establishing that the estimators developed in Section 4 are sound under the reasonable conditions. Before discussing our implementation and presenting the experimental results, we comment on the error probability in Theorem 9. This probability is small when each of the terms in the sum are small. Since each of these terms decays (at least) exponentially with $\delta$ times an expectation that grows linearly with the length of the string (assuming the mutation parameters are fixed), they will all generally be small.

Our estimators for the three rates require counting the number of $k$-spans with single deletion, $𝒟$, and the number of $k$-spans with no mutation, $𝒩$. Counting $𝒟$ and $𝒩$ can be challenging, particularly because $k$-mers do not contain the contextual information, and so we do not have access to their corresponding $k$-spans. An additional layer of complexity comes into play from the fact that identifying a $k$-mer with no mutation (or a single kind of mutation) is more difficult, considering many edge cases that may arise from inserting the same character that has been deleted. These challenges are cimcumvented when we use $k$-spans, and therefore, counting $𝒟$ and $𝒩$ solely from $k$-mers is not trivial and can be considered an interesting problem in and of itself. We therefore implemented an ad hoc solution to estimate $𝒟$ and $𝒩$ given the two strings $S$ and $S^{\prime }$. The steps for estimating $𝒟$ and $𝒩$ are as follows.

We start by extracting all $k$-mers in $S$, and building a de Bruijn graph using these $k$-mers using the cuttlefish tool [9]. We then extract the unitigs from this graph. Let the set of unitigs computed from $S$ be $U$. We also compute the unitigs in $S^{\prime }$ in a similar manner and call this $U^{\prime }$. We next take an arbitrary unitig $u$ from $U$, and align every unitig in $U^{\prime }$ with $u$. To allow for partial overlap, we use semi-global alignment by using the infix option in edlib [18], which makes sure that gaps at the beginning and at the end of the alignment are not penalized. For a particular $u\in U$, we align every $u^{\prime }\in U^{\prime }$ to make sure all relevant alignments are considered. We use these alignments to look at all windows of length $k$, and count $𝒟$ and $𝒩$ accordingly for $u$. We repeat this for all $u\in U$, and accumulate the measurements from individual $u$’s into a single global count.

The motivation behind using unitigs is that if there is an isolated mutation, and if the mutation is in the first or the last position of a $k$-mer, then there is no way to understand if the mutation is a substitution, an insertion, or a deletion only from the $k$-mers. The only way to resolve this ambiguity (and other similar ambiguities) is to scan beyond the context of $k$ characters – and unitigs are a natural way to do this. The core goal of our implementation of these steps described above was not to make it efficient, but rather to obtain a working solution. We found that executing these steps estimates $𝒟$ and $𝒩$ reasonably well, and the estimated rates are also acceptable. As such, we leave finding an efficient way to compute $𝒟$ and $𝒩$ as an open research question.

In this section, we present a series of experiments to evaluate the performance of the estimators detailed in Section 4. As discussed earlier, these estimators are sensitive to several input parameters, including $k$-mer size, sequence length, and the fraction of “A” characters in the sequence. Sections 7.1 through 7.3 explore the sensitivity of the estimators with respect to these parameters. In Sections 7.4 and 7.5, we estimate mutation rates across a wide range of known rate combinations. And finally, in Section 7.6, we demonstrate that our estimated substitution rate outperforms estimates obtained under a substitution-only mutation model. For the experiments in Sections 7.1 through 7.4, the original sequence is a randomly generated synthetic sequence. In these cases, we compute the number of $k$-spans containing a single deletion and the number of $k$-spans with no mutation directly from the known mutation process. For the experiments in Sections 7.5 and 7.6, we use real reference genomes as the original sequences. In these cases, the two types of $k$-span counts are estimated using the steps described in Section 6.

We have presented a mutation model that accommodates single-nucleotide substitutions, as well as insertions and deletions while retaining enough mathematical structure to admit closed-form rate estimators derived solely from $k$-mer statistics. From this model, we obtained algebraic estimators for the three elementary mutation rates: $p_s$, $p_d$, and $d$; we also proved a relatively tight sub-exponential concentration bound on $p_s$ that guarantees consistency as sequence length grows. We also identified regimes in which the estimation becomes ill-conditioned (i.e. large $k$, $p_d=0$, or sequence composition with $25\%$ “A”). These results establish a bridge between sequence evolution and combinatorial word statistics, thus providing additional tools for theoretical algorithmic computational biology.

In our prototype implementation, we demonstrated that our estimates remain accurate on simulated evolution of real genomes, and outperforms a substitution-only simple mutation model by avoiding spurious attribution of indel signals. While naive counting of unmutated and single-deletion $k$-mers sufficed to show practical accuracy of our estimators, this raises an interesting open problem: estimating the number of these unmutated and single-deletion $k$-mers efficiently for large scale data sets.

Several directions invite further investigation. First, incorporating the count $𝒮$ of single-substitution $k$-spans may illuminate why $p_s$ remains relatively stable even for moderately large $k$. Second, our framework can extend to heterogeneous or context-dependent rates by replacing global expectations with position-specific covariates. Third, coupling our estimators with sketch-based distance measures (such as in [8]) may provide a theory-backed avenue for larger scale applications such as phylogenetic placement in the presence of high indel activity. Finally, a more thorough investigation on real genomic data (where the unitig-based approach we used in the practical implementation starts to become infeasible) will be necessary to understand the utility of the mutation estimates in practice.

In summary, by utilizing probabilistic modeling and concentration inequalities, we provide a theoretical foundation and initial practical implementation for quantifying the parameters of a relatively complex mutation process directly from $k$-mers. We anticipate that these ideas will continue to inform new alignment-free computational biology tasks, particularly relevant as sequencing data continues to outpace traditional alignment-based paradigms.