A k-mer-Based Estimator of the Substitution Rate Between Repetitive Sequences

$K$-mer-based analysis of genomic data is ubiquitous. e.g. in genome assembly [1], error correction [2], read mapping [13], variant calling [29], genotyping [30, 7], database search [14, 9], metagenomic sequence comparison [26], and alignment-free sequence comparison [28, 20, 24]. One of the major challenges is the presence of repetitive $k$-mers, which adversely affects the practical performance as well as the theoretical analysis of downstream algorithms. One example is that heuristic read aligners like minimap2 [15] and even more rigorous ones like Eskemap [25] filter out highly repetitive $k$-mers in order to avoid explosive run times. Another example is the recent paper [27] that proved that sequence alignment can on average be done in almost $𝒪(n\log n)$ time but could not account for sequences with a high number of repeats.

One of the major advantage of $k$-mer-based methods is that they lend themselves more easily to sketching [16, 22], which is important for scaling to large-scale data. The groundbreaking Mash paper [20] was able to estimate the mutation rate between two genomes fast enough to be able to construct a phylogeny of 17 primate species in a tiny fraction of the time it would take an alignment-based method. Their approach uses an estimator based on the Jaccard similarity between the $k$-mer sketches of two sequences. However, the derivation behind their estimator assumes that the genomes have no repeats, making it inaccurate in highly repetitive regions. Other methods for estimating mutation rates are not designed for and/or not tested on highly repetitive sequences [34, 10, 18, 23].

In this paper we tackle the challenge of accounting for repeats when estimating the mutation rate. We assume that a string $t$ is generated from a string $s$ through a simple substitution process [5], where every nucleotide of $s$ mutates with a fixed probability $r$. Given the number of shared $k$-mers between $s$ and $t$ and the $k$-mer abundance histogram of $s$, we define our estimator $\widehat{r}$ as the solution to a polynomial equation, which can be solved using Newton’s method. We give a theorem to bound its bias, in terms of properties of $s$ (Theorem 3). Our estimator is designed to capture the most salient properties of the repeat structure of the genome, with the rest of the information being captured in the bias bounds. As a result, a user can decide a priori whether to trust our estimator, based on the quality of the bias bounds and on another heuristic we provide (Theorem 4).

We evaluate our estimator $\widehat{r}$ empirically on various sequences, including the alpha satellite centromeric region of human chr21 and the highly repetitive human RBMY gene. For such repetitive sequences, our estimator remains highly accurate, while the repeat-oblivious estimator of the kind used by Mash is unreliable. We make a comprehensive evaluation of $\widehat{r}$ across the spectrum of $k$ and $r$ values, which can guide a user towards choosing a $k$ value for their analysis. We also show that our estimator can be used on top of a FracMinHash sketch, without systematically effecting the bias. Our software is available on GitHub [32].

Let $s$ be a string and let $k>0$ be a parameter indicating the $k$-mer size. We will index string positions from 1. We further assume in this paper that $s$ is at least $k$ nucleotides long. We use $L$ to denote the number of nucleotides in the string minus $(k-1)$, or, equivalently, the number of $k$-mers in $s$. For $1\le i\le L$, let $s_i$ be the $k$-mer starting at position $i$ of $s$. Let $s{p}^k(s)$ be the set of all distinct $k$-mers in $s$, also called the $k$-spectrum of $s$. We let $L_0$ be the size of $s{p}^k(s)$, i.e. the number of distinct $k$-mers in $s$. Given a $k$-mer $\tau$, we will use the shorthand $\tau \in s$ to mean $\tau \in s{p}^k(s)$. Given two strings $s$ and $t$, we define $I(s,t)\triangleq |s{p}^k(s)\cap s{p}^k(t)|$ as the number of $k$-mers shared between them. We will usually use the shorthand of $I$ for $I(s,t)$. Given two $k$-mers $\tau$ and $\nu$, we use $\text{HD}(\tau ,\nu )$ to denote their Hamming distance.

Let $K$ be a set of $k$-mers and let $s$ be a string. We let $occ(K)$ denote the number of positions $i$ in $s$ such that $s_i\in K$. When $K$ consists of a single element $\tau$, we simply write $occ(\tau )$. A set of positions $𝒥$ is said to be a set of occurrences of $K$ if for all $i\in 𝒥$, we have $s_i\in K$. A set of occurrences is said to be non-overlapping if, for all distinct $i,j\in 𝒥$, $|j-i|\ge k$. We let $sep(K)$ be the maximum size of a set of non-overlapping occurrences of $K$, also referred to as the separated occurrence count. Observe that $0\le sep(K)\le occ(K)$. The abundance histogram of a string $s$ is the sequence $(a_1,\mathrm{\dots },a_L)$ where $a_i$ is the number of $k$-mers in $s{p}^k(s)$ that occur $i$ times in $s$. Note that $L_0={\sum }_{i=1}^L{a}_i$.

We will consider the following random substitution process, parameterized by a rate $0\le r\le 1$. Given a string $s$, it generates an equal-length string where, independently, the character at each position is unchanged from $s$ with probability $1-r$ and changed to one of the three other nucleotides with probability $r/3$.

In this paper, we address the following problem. Let $0\le r\le 1$ be a substitution rate. Let $s$ be a string and let $t$ be generated from $s$ using the substitution process parametrized by $r$. Let $I_{\text{obs}}=I(s,t)$ be the observed spectrum intersection size. Given $I_{\text{obs}}$ and the abundance histogram of $s$, the problem is to estimate the mutation rate $r$.

The bias of an estimator $\widehat{r}$ for $r$ is defined as $\mathrm{E}[\widehat{r}]-r$. A good estimator should have a small absolute bias, one that is within the error tolerance of downstream applications. For our problem, directly finding an estimator for $r$ with provably small bias turned out to be technically challenging. Instead, we provide an estimator for $q\triangleq 1-{(1-r)}^k$, which corresponds to the probability that a $k$-mer occurrence contains at least one substitution. There is a natural one-to-one correspondence between an estimator $\widehat{q}$ of $q$ and an estimator $\widehat{r}$ of $r$ via the equation $\widehat{q}=1-{(1-\widehat{r})}^k$. Thus, an alternative to bounding the bias of $\widehat{r}$ is to bound that of $\widehat{q}$; i.e., bound $\mathrm{E}[\widehat{q}]-q=\mathrm{E}[1-{(1-\widehat{r})}^k]-(1-{(1-r)}^k)$. While the difference between the two approaches may intuitively seem minor, bounding the bias of $\widehat{q}$ turned out to be more technically feasible.

The only previously known estimator for this problem is what we refer to as the repeat-oblivious estimator¹¹1 What we describe is based on the estimators used in [8, 20, 12], but with two important differences. The first is that we use the modification adopted in the follow up work of [24] and described in Appendix A.6 of [3]. The second is that the original estimator was calculated from the Jaccard similarity between two sequences; however, under our substitution process model, we can state more simply in terms of the spectrum intersection size. and denote by $r_{\text{obl}}$. The derivation of this estimator assumes that 1) $s$ has no repeats and 2) if a $k$-mer mutates, it never mutates to anything that is already in $s$. The estimator is then derived using a simple two step technique, called the method of moments [6]. The first step is to derive $\mathrm{E}[I]$ in terms of $r$. Under the assumptions (1) and (2), $\mathrm{E}[I]=L{(1-r)}^k$. The second step is to take the observed value $I_{\text{obs}}$, plug it in place of $\mathrm{E}[I]$, and solve for $r$. In this case, one solves the equation $I_{\text{obs}}=L{(1-r)}^k$ and gets the estimator $r_{\text{obl}}=1-{\left(\frac{I_{\text{obs}}}{L}\right)}^{1/k}$ and the corresponding estimator $q_{\text{obl}}=1-\frac{I_{\text{obs}}}{L}$. The estimator $q_{\text{obl}}$ is an unbiased estimator of $q$ when (1) and (2) hold, but the bias for general sequences is not known. We are also not aware of any results about the bias of $r_{\text{obl}}$, even under assumptions (1) and (2). As we will show in Section 6, the repeat-oblivious estimator has a large empirical bias when the assumptions are substantially violated.

On the one hand, the repeat-oblivious estimator does not at all account for the repeat structure of $s$. On the other hand, an estimator that would fully account for it seems to be challenging to derive, analyze, and compute. Moreover, such an estimator would likely be superfluous for real data. Instead, our approach is intended to achieve a middle ground between accuracy and complexity by accounting for the most essential part of the repeat structure in the estimator and expressing the non-captured structure in the bias formula. We will show that under assumption (2) and the assumption that all $k$-mer occurrences are non-overlapping in $s$,

Following the method of moments approach, we define our estimator $\widehat{r}$ as the unique solution (the uniqueness is shown in Lemma A.2) to this equation when plugging in $I_{\text{obs}}$ in place of $\mathrm{E}[I]$. Though we are not able to analytically solve for $\widehat{r}$, we can find the solution numerically using Newton’s Method.

Note that the assumptions we make are not necessary to compute $\widehat{r}$ and only represent the ideal condition for our estimator. Our theoretical and experimental results will quantify more precisely how the deviation from our assumptions is reflected in the bias.

Recall that we define $\widehat{q}=1-{(1-\widehat{r})}^k$ and, as mentioned earlier, we will prove the theoretical results on the bias of $\widehat{q}$, rather than $\widehat{r}$. First, we need to derive the expectation and variance of the intersection size. A closed-form expression for even the expectation is elusive, so we will instead use an approximation and derive bounds on the error. The idea behind our bounds is that the error becomes small on the types of sequences that occur in biological data.

We want to underscore that when we make probabilistic statements, it is with regard to the substitution process and not with regard to $s$. We do not make any assumptions about $s$, and, in particular, we are not considering the situation where $s$ itself is generated randomly.

First, it is useful to express $I\triangleq I(s,t)$ as a sum of indicator random variables. Let us define $E_i^{\tau }$ as event that $t_i=\tau$ and $E^{\tau }={\cup }_{i=1}^L{E}_i^{\tau }$ as the event that at least one position in $t$ contains $\tau$. By linearity of expectation, we have

Let $ℱ(q)\triangleq L_0-{\sum }_{i=1}^L{a}_i{q}^i$. Relying on the approximation $\mathrm{E}[I]\approx ℱ(q)$ (i.e. Equation 1), we define $\widehat{q}$ as the solution to $I_{\text{obs}}=ℱ(q)$, or, equivalently, $\widehat{q}={ℱ}^{-1}(I_{\text{obs}})$. We show that this approximation holds when we assume that 1) $\mathrm{Pr}[E_i^{\tau }]=0$ when $s_i\ne \tau$ (see footnote²²2 We note that this assumption is not theoretically precise, because forbidding a $k$-mer at position $i$ from mutating to $\tau \in s$ usually implies that there is at least one $\nu \notin s$ that the $k$-mer at position $i+1$ can no longer mutate to. Because of these dependencies, there are downstream effects on the probability space that are complex to track. A theoretically robust alternative was given in [5] via the $k$-span formulation of the problem. It could be used to formalize the assumption here, however, in this paper, we only use the assumption to give intuition for the estimator and do not use it in any formal theorem statements or proofs. ) and 2) all occurrences of $\tau$ are non-overlapping in $s$:

The underlying philosophy for our estimator is that while these assumptions are not perfectly satisfied on real data, in most cases the contribution due to violations of these assumptions is small. To make this mathematically precise, we will bound the difference between $\mathrm{E}[I]$ and $ℱ(q)$ in terms of an expression that can be calculated for any $s$.

The difference between $ℱ(q)$ and $L_E$ (i.e. ${\sum }_{\tau \in s}{q}^{sep(\tau )}-q^{occ(\tau )}$) is close to 0 when the number of $k$-mers with overlapping occurrences is close to 0. On the other hand, the difference between $ℱ(q)$ and $U_E$ (i.e. ${\sum }_{\tau \in s}{\beta }_{\tau }$) is never zero (except in corner cases). However, the largest terms contributing to this difference are due to pairs of non-identical $k$-mers that have a small Hamming distance to each other. Thus, the difference becomes small when the number of “near-repeats” is small.

Next, we upper bound the variance.

Theorem 2 gives an upper bound on $\mathrm{Var}[I]$ in two forms. The first one is more precise because it is a function of $\mathrm{E}[I]$. However, since we are not able to compute $\mathrm{E}[I]$ exactly, the second form allows us to plug in the upper bound on $\mathrm{E}[I]$ from Theorem 1.

Given the bounds on $\mathrm{E}[I]$ and $\mathrm{Var}[I]$, we are now able to bound the bias of $\widehat{q}$.

The derivatives of $ℱ(q)$ have straightforward closed-form expressions, since $ℱ$ is a polynomial in $q$. We do not have a closed-form solution for $f$, but it can be evaluated numerically using Newton’s method. Thus, for any given sequence $s$, we can precompute the bounds of our $\widehat{q}$ estimator bias for any value of $q$. Due to space limitations, we do not further elaborate on the algorithm to compute the bounds in Theorem 3 or on its runtime analysis.

When the observed intersection is empty, there is a loss of signal and it becomes challenging for any intersection-based estimator to differentiate the true substitution rate from 100%. The following theorem gives an upper bound on the probability that the intersection is empty, as a function of $L,k$, and $r$. In Section 6, we will show how it can be used to make a conservative decision that the computed estimate is unreliable.

This section contains the proofs of our theoretical results. In particular, we will prove Theorems 1, 2, and 3 from the previous section. The proof of Theorem 4 is left for the Appendix. We start by proving a couple of preliminary facts that will be used in the proofs of these theorems. First, we consider the probability of the event $E_i^{\tau }$, which is straightforward to derive.

In this section, we evaluate the empirical accuracy and robustness of our estimator, used by itself or in combination with sketching.

In this paper, we propose an estimator for the substitution rate between two sequences that is robust in highly repetitive regions such as centromeres. Our experiments validated its performance across a broad range of $k$ and $r$ values. We provide theoretical bounds on our estimator’s bias (specifically on the bias of $\widehat{q}$), and show that it accurately captures the estimator’s empirical mean in most scenarios.

For large values of $k$ and $r$, i.e., when $q$ is large, the intersection of the $k$-spectra tends to be empty with high probability, which is a foreseeable limitation for all intersection-based estimators. To address this, we introduce a heuristic criterion, $P_{empty}$, which depends only on the number of $k$-mers $L$, the $k$-mer size $k$, and the substitution rate $r$. This criterion allows us to heuristically identify parameter settings under which the estimator becomes unstable.

We also showed how our estimator can be easily combined with FracMinHashing. Empirical results show that sketching does not introduce systematic bias, albeit at the cost of increased variance.

We do not perform a runtime analysis of our estimator because it completes in less than a second on our data. The runtime of our estimator is the time it takes to solve an equation numerically using Newton’s method. Since $ℱ(q)$ is a polynomial and the solution is constrained to the interval $[0,1]$, Newton’s method converges in $𝒪(\log \log (1/ϵ))$ iterations, where $ϵ$ is the target precision. Each iteration involves evaluating $ℱ$ and its derivative, which takes time proportional to the number of non-zero $a_i$ terms. Except for esoteric corner cases, the number of such terms is small in practice.

The immediate open problem is to tighten the theoretical bounds on the bias. Future work could thus focus on deriving a tighter variance bound to strengthen the theoretical characterization of $\mathrm{E}[\widehat{q}]$. A bigger open question is how to derive confidence intervals. This is a more challenging problem than bounding the bias because it requires a deeper understanding of the estimator’s distribution.

Our estimator could potentially be extended to work on unassembled sequencing reads, as opposed to assembled genomes. Our method does not rely on the $k$-mer multiplicities in the intersection size, making it amenable to such a scenario. Still, one of the limitations of our estimator is the need to know the abundance histogram of the source string. A tool like GenomeScope [31] can estimate the abundance histogram from sequence data $k$-mer counts. Alternatively, the user may choose to use an abundance histogram from a related genome, as related genomes are likely to have similar abundance histograms. Fully adapting this estimator to work with sequencing data remains an important future work.