Proving sequence aligners can guarantee accuracy in almost O(m log n) time through an average-case analysis of the seed-chain-extend heuristic

Our contribution

Our goal in this paper is to close the gap between theory and practice, rigorously justifying the heuristics used in some of the most widely used alignment software. To this end, we turn to the methods of average-case analysis (Szpankowski 2001), which give us a way of breaking through the quadratic barrier of alignment (Medvedev 2022b). Given a random string, we define a substitution model giving rise to a distribution on pairs of inputs and then average our analysis on pairs of strings over this distribution.

Recently, Ganesh and Sy (2020) also used probabilistic analysis to show that a heuristic algorithm based on banded alignment can run in O(n log n) time for two length-n sequences and is optimal with high probability. However, their specific method has not yet seen any usage in practical software, and their analysis is invalid in the case of read mapping as it only pertains to two nearly-equal-length strings. Thus, we turn instead to the analysis of an empirically battle-tested heuristic for sequence alignment: seed-chain-extend. Seed-chain-extend is a well-established technique for comparative genomics (Myers and Miller 1995; Abouelhoda and Ohlebusch 2005), and recently, the addition of sketching (or subsampling) has made it popular for long-read aligners (Chaisson and Tesler 2012; Sović et al. 2016; Ren and Chaisson 2021), including minimap2 (Li 2018), the primary algorithm our model of seed-chain-extend is based on.

We provide, to the best of our knowledge, the first average-case bounds on runtime and optimality for the sketched k-mer seed-chain-extend alignment heuristic under a pairwise mutation model. Our optimality result shows that for large-enough k-mer size k, the alignment is mostly constrained to be near the correct diagonal of the alignment matrix, and that runtime is close to linear when the mutation or error rate is reasonably small. We also show that subsampling $\mathrm{\Theta }\left({\displaystyle \frac{1}{\log n}}\right)$ of k-mers asymptotically reduces our bounds on chaining time but not extension time. Our results give a theoretical justification for both the empirical accuracy and subquadratic runtime of seed-chain-extend.

Mutation model

Let S = x₁x₂…x_n+k−1 be a random uniform string with n + k − 1 i.i.d. letters on an alphabet of size σ for some k = Θ (log n). Let S′ = y_p+1y_p+2…y_p+m+k−1 be a substring of S of length m + k − 1 starting at a fixed position p + 1 with each character independently mutated to a different letter with probability θ. Although notationally a little confusing at first, the k − 1 term ensures S, S′ contain exactly n and m k-mers, respectively: k-mers are in many ways the natural unit of measurement rather than individual characters. We model only point substitutions here and not indels. Independent substitution models have been considered in theoretical work (Ganesh and Sy 2020; Blanca et al. 2022; Shaw and Yu 2022) and, importantly, also have shown to be useful empirically (Chaisson and Tesler 2012; Ondov et al. 2016; Sarmashghi et al. 2019). Also, although genomes can be repetitive, on the level of k-mers, a random model has been shown to be reasonable (Fofanov et al. 2004). We discuss other possible random models in the Discussion section.

Modeling seed-chain-extend

A brief overview of seed-chain-extend-based alignment is given as follows: First, a subset of k-mers in both S, S′ is taken as seeds, and exact seed matches between S, S′, called anchors, are obtained. We only use k-mer seeds in this study, although other types of seeds are possible (Keich et al. 2004). An optimal increasing subsequence of possibly overlapping anchors based on some score is then collected into a chain, in which increasing is defined with the standard precedence relationship (Jain et al. 2022) between k-mer anchors (see subsection “Chaining” below) (see Fig. 5A, below). The chain is extended into a full alignment by aligning between anchor gaps in the chain.

Extension and chaining runtimes

Given sorted anchors, let T_Chain be the time spent finding an optimal chain. T_Chain depends on the objective function (Abouelhoda and Ohlebusch 2005; Otto et al. 2011; Mäkinen and Sahlin 2020; Jain et al. 2022). Because our gap costs are linear, T_Chain = O(Nlog N), where N is the number of anchors (Abouelhoda and Ohlebusch 2005). For extension time T_Ext, let $\left(G_1,{G^{\mathrm{\prime }}}_1\right),\dots ,\left(G_{u-1},{G^{\mathrm{\prime }}}_{u-1}\right)$ be the size of the gaps for an optimal chain. G_ℓ indicates the length of the substring between the k-mers i_ℓ, i_ℓ+1 on S and $G_{\ell }^{\mathrm{\prime }}$ similarly for S′; $G_{\ell },G_{\ell }^{\mathrm{\prime }}$ can be zero. The extension time is $T_{Ext}={\sum }_{\ell =1}^{u-1}O\left(G_{\ell }{G^{\mathrm{\prime }}}_{\ell }\right)$. Using the fact that $\sum G_{\ell }\le n$ and $G_{\ell }^{\mathrm{\prime }}\le \sum G_{\ell }^{\mathrm{\prime }}\le m$, one can get that $T_{Ext}={\sum }_{\ell =1}^{u-1}O\left(G_{\ell }m\right)=O\left(nm\right)$. We will show that the expected runtime is better than this upper bound, but it serves as a useful worst case. Because S, S′ are random strings, T_Chain, T_Ext, N, and the alignment itself all become random variables. Our goal will be to bound E[T_Chain] and E[T_Ext].

Theoretical results

First, a few definitions are in order. Recall that $\left|S\right|=n+k-1\sim n$, and $\left|S^{\mathrm{\prime }}\right|=m+k-1\sim m$ are our string lengths. We define σ > 1 to be the size of the alphabet and 0 < θ < 1 to be the probability a base mutates. Our theory holds for any alphabet size, but we use σ = 4 for specifying numerical constants. Let $\log =\mathrm{lo}{\mathrm{g}}_{\sigma }$ with base σ, and let ln be the natural logarithm. Let k = C log n for a fixed C > 0 so that key quantities can be expressed interchangeably as σ^k = n^C and (1 − θ)^k = n^−Cα, where $\alpha =-{\log }_{\sigma }\left(1-\theta \right)>0$. We define the actual goodness of the chain in terms of recoverability in the Methods subsection “Homology and recoverability”; it measures the fraction of homologous bases under our mutation model that could potentially be recovered by extending through an optimal chain and only depends on the chain, not the actual alignment.

The above condition on C is equivalent to $C\alpha ={\displaystyle \frac{2\alpha }{1-2\alpha }+\delta }$ for any δ > 0. In practice, we will make δ very small so k is not too large. ${\displaystyle \frac{2\alpha }{1-2\alpha }}$ is plotted in Supplemental Figure S7; it is convex and thus leads to the bound Cα < 2.43 · θ. For example, if |S| = |S′| and θ = 0.05, which is approximately the level of divergence seen between human and chimpanzee genomes (Britten 2002), ${\displaystyle \frac{-2{\log }_4\left(0.95\right)}{1+2{\log }_4\left(0.95\right)}<0.08}$, so for Cα = 0.08, the running time is O(n^1.08log (n)). Even for relatively large values of n, say the entire human genome, which has size of approximately 3 billion, n^1.08 is essentially linear as (3,000,000,000)^0.08 ≈ 5.73.

Our proof follows in three main steps. First, we bound the first and second moments of the random variable N, which denotes the number of anchors, implying that chaining is fast. Second, we use concentration inequalities for sums of dependent random variables and exploit the structure of chaining to show that with high probability, an optimal chain does not deviate much from the chain with only “homologous anchors” (Fig. 5A). The failure probabilities will be on the order of $\mathrm{\Theta }\left({\displaystyle \frac{1}{n}}\right)$, allowing us to also bound everything in expectation. Lastly, we bound the expected runtime of extension through gaps between homologous anchors.

Asymptotically, the chaining runtime is smaller than the extension runtime. However, the implied constants can be much smaller for the extension term (Kalikar et al. 2022), so it is practically useful to reduce the runtime of chaining via sketching. Let 0 < 1/c ≤ 1 be the expected fraction of selected k-mers for a sketching method. We use the open syncmer k-mer seeding method (Edgar 2021), which has a useful mathematical property (Theorem 7), giving the following:

This shows that the asymptotic upper bound on extension runtime is the same even if we let c grow with n like c = Θ (log n) < k, leading to the following conclusion: Sketching can reduce chaining time without increasing extension time much. Other seeding schemes used, for example, minimizers (Li 2018; Rautiainen and Marschall 2020; Colquhoun et al. 2021; Sirén et al. 2021) or FracMinHash (Irber et al. 2022), behave differently, but our techniques provide intuition and, in some cases, can be extended.

Simulated genome alignment experiments

Although the model of seed-chain-extend we consider in the theory is based off of real aligners (e.g., minimap2), real aligners implement many additional tricks to make things work better. As such, to empirically validate our theory, we implemented a basic version of sketched seed-chain-extend and tested it on simulated random sequences with independent point substitutions. For the chaining step, we implemented an AVL tree–based max-range-query method (Mäkinen et al. 2015; Li et al. 2020). For the extension step, we used a standard dynamic programming (DP) algorithm implemented in rust-bio (Köster 2016).

We let k = C log n, where $C={\displaystyle \frac{2}{1-2\alpha }}$, and do seed-chain-extend on two length n sequences in this section, but we do an additional experiment for a substrings of length m = n^0.7 + 100 in Supplemental Figure S10. As k cannot be fractional, we let $n_k=4^{\frac{k}{C}}$ for varying integer values of k. We will set c = k − 7 = C log (n) − 7, which grows with n. The constant 7 is arbitrary but chosen to be sufficiently large so that s-mers in the open syncmer method are not so small as to make the method degenerate. We found that recoverability was always quite high and that breaks did not occur very often in actual simulations; we show this in Supplemental Figure S9. Thus, the primary focus of the empirical results will be runtime.

Accuracy of asymptotic extension runtime predictions

We first empirically investigate our upper bound on expected extension runtime, which is O(n^1+Cαlog (n)) for both sketched and nonsketched extension when c = Θ (log n). Assuming that the runtime is simply $\lambda n_k^{1+C\alpha }\log \left(n_k\right)$ for some fixed constant λ, we can predict the ratio of the runtimes as ${\displaystyle \frac{n_{k+1}^{1+C\alpha }\log \left(n_{k+1}\right)}{n_k^{1+C\alpha }\log \left(n_k\right)}}$. Of course, this is incorrect for small n, as smaller terms may dominate runtime, but we expect it to be reasonably accurate for large n. We plot the empirical and predicted ratios of extension runtimes in Figure 1A.

Figure 1.

Runtime results on seed-chain-extend between for two sequences of length n and $k=C\log n={\displaystyle \frac{2}{1-2\alpha }\log \left(n\right)=9,10,\dots ,19}$, where θ = 0.10 and α = −log ₄(1 − θ). (A) The y-axis shows the runtime ratios (with 95% confidence intervals) for iteration k + 1 divided by runtime for iteration k, where the sequence length n_k is plotted on the x-axis. As n grows, both the sketched and nonsketched extension ratios asymptotically approach the predicted ratio of $n^{1+\frac{2}{1-2\alpha }\alpha }\log \left(n\right)$ runtimes. (B) Multiplicative speed-up for sketched versus nonsketched alignment when sketching with density ${\displaystyle \frac{1}{c}={\displaystyle \frac{1}{k-7}}}$, where k = C log n. The slow-down for extension flattens out, whereas the chaining speed-up is almost linear on the log-scaled x-axis; chaining speed-up dominates extension slow-down as predicted by our theory.

Open in new tabDownload PowerPointLink to figure

Figure 1A shows that our upper bound looks reasonable for both sketched and nonsketched cases for θ = 0.10. In Supplemental Figure S8, we show the same plots for θ = 0.05. The empirical results never cross the predicted ratio, but this is not too surprising as the predicted ratio is only approximate. Importantly, the empirical extension runtime ratios slope downward, which agrees with the prediction.

Real nanopore read alignment experiments

To test the applicability and generalizability of our theorems to real, nonsimulated data, we performed an experiment by aligning real long-read data from Oxford Nanopore Technologies (ONT) using seed-chain-extend. We use the sketched seed-chain-extend aligner in the simulated experiments with a slight modification, which is explained in the next section. Importantly, our two main assumptions with respect to uniform random strings and point substitutions are violated in this setup: Real genomes are not uniformly random strings, and ONT reads contain a significant amount of indel errors (Delahaye and Nicolas 2021).

Our data set consists of five reference genomes and corresponding long-read data sets. The references consist of a viral genome (SARS-CoV-2), a bacterial genome (Escherichia coli), a fungi genome (Magnaporthe oryzae), an insect genome (Drosophila melanogaster), and a human genome (Homo sapiens). For each genome, we aligned a corresponding set of publicly available ONT long-reads, which may have different length and error distributions. The data sets can be found in Supplemental Table S1.

We are interested in testing the predicted results as a function of m (the read length) and n (the genome size). We let $k=C\log n={\displaystyle \frac{2}{1-2\alpha }\log n}$ as before. Because α and C in our results depend on the error rate θ, we first aligned each read set to the genomes with minimap2 and retained only reads with gap-compressed identity of > 94.5% and < 95.5% so that θ ≈ 0.05 is fixed. In addition, we only used reads for which minimap2 aligned >90% of the read, to avoid benchmarking on unalignable reads.

Evaluating model assumptions

We first introduced the following filters to our basic seed-chain-extend alignment algorithm for computational reasons:

1. We filtered reads when the number of anchors is more than 10 times the number of bases. This occurs because of extremely repetitive k-mers and can cause the chaining step to stall.

2. We did not consider reads when the chains had gaps of a length >10,000 bp, which could occur because of structural variations or failure of chaining and can cause extension to stall.

The two filters also measure how significantly our model assumptions, with respect to independent substitutions and uniform random strings, are violated. We report the number of reads that fail the above filters in Table 1. Only a small fraction of reads fails the 10,000-bp gap filter, but the repetitive k-mer filter is violated more frequently as the genomes become more complex. Still, >80% of the human reads and > 97% of the D. melanogaster reads are not deemed too repetitive. Mapping to repetitive regions is an active area of research (Jain et al. 2020), which is beyond the scope of this paper, but effective heuristics such as repetitive k-mer masking (Li 2018; Sahlin 2022) exist.

Table 1.

Counting the number of reads that did not pass our gap filter and repetitive k-mer filter

Open in new tabDownload PowerPointLink to table

Extension and chaining runtimes are well predicted by theory

In Figure 2, we plot the times of chaining and extension and perform a robust linear regression using the Siegel estimator (Siegel 1982) as a function of the read length m. For this plot specifically, we only timed alignments in which the read was > 90% aligned by our seed-chain-extend implementation. For each plot, the genome size n is held fixed. For fixed n, all runtimes are well approximated by a linear function in m as predicted by our theory, although chaining times for the human genome have higher variance owing to the presence of repetitive k-mers.

Figure 2.

Extension and chaining time for real ONT reads with ≈5% sequence divergence. k-mer size was chosen to be k = C log n for reference length n and $C={\displaystyle \frac{2}{1-2\alpha }\approx 2.16}$ with error parameter θ = 0.05. Exact data sets are described in Supplemental Table S1. Note the scaled axes for the SARS-CoV-2 reads, which were much smaller than the other data sets. The well-fit linear regression lines indicate essentially linear runtime in read length with fixed reference length n (i.e., fixed k = C log n and constant C), although larger human chaining time variance is owing to repetitive k-mers.

Open in new tabDownload PowerPointLink to figure

In Figure 3, we plot the slopes of the regressions for extension time obtained from the first row of Figure 2, which measures the dependence on n given fixed m. This dependence is reasonably fit by a log n function (R² = 0.766), so O(mlog n) is not a bad approximation for the runtime in practice. However, our theory leads to a n^Cαlog n dependence, which comes from the big O extension runtime. Cα ≈ 0.08 when θ = 0.05 in our setup, so by using a n^0.08log n function instead, we end up getting a better fit (R² = 0.928), as predicted by our theory.

Figure 3.

Data points represent the slopes of the linear regressions in Figure 2 for extension, with the corresponding value of k (which is k = C log n for a constant $C={\displaystyle \frac{2}{1-2\alpha }}$). This dependence of the genome size, n (x-axis), is decently approximated by a naive A₁log (n) + B₁ fit, where A₁ and B₁ are parameters. However, our theory states that the dependence should be log (n)n^Cα with Cα ≈ 0.08 when θ = 0.05. Fitting A₂log (n)n^0.08 + B₂ gives a better R² value (0.928 vs. 0.766) with the same number of parameters (two parameters for both fits), indicating the goodness of our theoretical predictions.

Open in new tabDownload PowerPointLink to figure

Recoverability upper bound by aligned fraction

Recoverability is not measurable on real, nonsimulated reads because true sequence homology is not known. However, we can upper bound recoverability by simply measuring the length of the chain relative to the read length, as any homologous bases outside of the chain (i.e., not within the first and last anchor) cannot be recovered under our model. We call this the aligned fraction, which our theory also predicts should be $>1-O\left({\displaystyle \frac{1}{\sqrt{m}}}\right)$ in expectation (with implied constants depending on θ and n). Practically, the aligned fraction is easily interpretable, and the ability to predict the aligned fraction is useful for ensuring that reads with a given error rate and k-mer size are long enough to align to a reference.

We also observed that the aligned fraction was the main cause of recoverability loss in the simulated experiments (see Supplemental Fig. S9), and it turns out that the main $O\left({\displaystyle \frac{1}{\sqrt{m}}}\right)$ term in the theoretical recoverability result exactly comes from bounding the aligned fraction in the proof (Supplemental Lemma S7), so this quantity is theoretically relevant as well.

In Figure 4, we plot the aligned fraction for each data set. We perform a least-squares fit of the function $1-{\displaystyle \frac{a}{m^b}}$, where a and b are parameters. We exclude points with < 0.25 aligned fraction if the read length is greater than 1000 when fitting because these points appear to be outliers caused by poor chaining and can bias the least-squares fit. The resulting fits for all plots are reasonable, with R² > 0.5 except for the D. melanogaster data set. The seemingly small value of R² can be explained by the fact that most reads (passing our filters) are well aligned, so the data are almost constant and ≈1. Constant data will have R² = 0, so we cannot do much better in this case. Nevertheless, b > 0.85 on all data sets, suggesting that our $\sqrt{m}$ term may be too conservative. However, the exact parameter values in the fit are highly sensitive to the fitting methodology, so the exact values displayed are only suggestive and not to be overinterpreted.

Figure 4.

The aligned fraction of the aligned reads, which is chain length divided by read length. We fit a function of the form $1-{\displaystyle \frac{a}{m^b}}$, where m is the read length (x-axis), and display the resulting R² values. Aligned fraction is an upper bound for recoverability, and our lower bound on recoverability (and also aligned fraction) is $1-O\left({\displaystyle \frac{1}{\sqrt{m}}}\right)$, which the data suggest may be too conservative.

Open in new tabDownload PowerPointLink to figure

Improving bounds for the current random model

For runtime, it seems unlikely that the expected runtime is truly quasilinear owing to the fundamental quantity (1 − θ)^k = (1 − θ)^{C log n} = n^−Cα, the likelihood of a k-mer match, decreasing faster than logarithmically in n when k = Θ(log n); this turns out to be the cause of the n^Cα term in the runtime (see Methods). Despite this, we believe there can be significant improvements for our bounds in Theorems 1 and 2. The most significant aspect that we believe can be improved is the restriction on the constant C. C is required to be $>{\displaystyle \frac{2}{1-2\alpha }}$, which seems unsatisfactory because this leads to high k-mer sizes, for example, k = 34 for θ = 0.05 on the human genome. This is much larger than the k-mer sizes used by noisy long-read aligners in practice (minimap2 defaults to k = 15).

The restricted values of C is because of our analysis of spurious anchors relying on a relatively weak variance bound (Lemmas 3 and 4). To use the bound effectively, C must be quite large. The variance bound is also the cause for the $1-O\left({\displaystyle \frac{1}{\sqrt{m}}}\right)$ bound for recoverability, which we also believe can be tightened based on both the simulated and real experiments. To surpass these bounds, a deeper understanding of spurious anchor subsequences must be developed. Spurious subsequences of anchors is a very similar problem to common subsequences in random strings for extremely large alphabets, a topic that has had much theoretical attention (Chvátal and Sankoff 1975; Navarro 2001; Kiwi et al. 2005; Lember and Matzinger 2009).

Generalizing the random model

Our current random model does not model complex genomes properly. It is well known that genomes can be extremely repetitive (de Koning et al. 2011), which violates the uniform random string assumption. This is why the expected runtime of chaining deviates from our predictions in Figure 2 for the human genome. One idea for modeling repeats is to use a k-mer Markov model for the random strings instead of an independent uniform string. However, the analysis of such models is highly nontrivial (Reinert et al. 2000), and our techniques for analyzing k-mers in this work may not apply. Regardless of how repeats are modeled, any sort of nontrivial analysis would still require a better handling of spurious anchor subsequences, which is exactly the issue that needs to be tackled in order to improve the bounds for our current model.

Another direction is to incorporate indels into the random model. The technical reason we opted for modeling point substitutions instead of indels is that indels complicate the statistics of k-mer matching. For example, an indel of the nucleotide A may not “change” the 4-mers in the sequence AAAAAA, and sequence homology becomes ambiguous. A formulation of an “indel channel” as a model of sequence mutation was used by Ganesh and Sy (2020) for analysis of edit distance, but the resulting analysis for this model was more complex than in the independent substitution case.

Generalizing our model to small indels seems less pressing as our theory was successful for real nanopore data, which have many small indel errors, showing that the substitution model already generalizes well to small indels. However, modeling large structural variations could be a relevant problem.

Homology and recoverability

Let the interval [a..b] be the discrete interval {a, a + 1, …, b} and S[a..b] be the substring indexed by [a..b]. S′ is a mutated version of the substring S[p + 1..p + m] of length m starting at p + 1, so the optimal alignment should map S′ to the homologous indices [p + 1..p + m]. We now define recoverability as the number of homologous bases one can possibly recover by seed-chain-extend; a visual example is shown in Figure 5.

We define $Align\left(\mathcal{C}\right)$ carefully in Section B in the Supplemental Materials and give a visual explanation in Figure 5B. Note that an optimal chain with respect to our linear-gap cost objective u − ζ[(i_u − i₁) + (j_u − j₁)] does not directly optimize for recoverability. We provably find the chain with the optimal linear-gap cost, and then we will argue that it leads to good recoverability.

Figure 5.

Seed-chain-extend visualized on an alignment matrix. Anchors are short diagonal matches in the matrix. Extension corresponds to performing dynamic programming (DP) on the submatrix between anchors. (A) k-mer anchors under mutations and their corresponding alignment matrix. Blue anchors are homologous anchors, and red are spurious anchors. (B) $Align\left(\mathcal{C}\right)$ consists of anchors and DP matrices. Recoverable bases correspond to the bases on the homologous diagonal where a k-mer lies or is accessible by the green DP matrix between gaps. Breaks cover nonrecoverable sections.

Open in new tabDownload PowerPointLink to figure

Traditionally, alignment optimality is defined using a generalized edit distance (Vingron and Waterman 1994). However, these distances are only used as a proxy for detecting sequence homology (Thorne et al. 1991; Batzoglou 2005; Lunter et al. 2005). We know the true underlying sequence ancestry in our model, so defining optimality with respect to sequence homology suits the actual goal of sequence alignment. The reason for the name “recoverability” is that extension could potentially recover all recoverable bases, but this depends on the extension algorithm and is not guaranteed (Lunter et al. 2008).

Under our model, the trivial O(1) alignment that aligns S′ back to the originating substring of S without indels is the most homologous. This may seem to make our results superfluous; however, remember that the algorithm does not know where S′ “begins” if m < n. Also, although we do not attempt the case with indels, seed-chain-extend is still valid when indels are present, whereas the trivial solution does not allow for indels.

We will lower bound $\mathbb{E}\left[R\left(\mathcal{C}\right)\right]$. To do this, we will work with a more natural object called a break.

We prove Lemma 1 in the Supplemental Materials. The concept of a break is illustrated in Figure 5B. Breaks cover nonrecoverable regions, so subtracting the breaks from the span of the anchors lower bounds the recoverability.

Fundamental tools and bounds

In this section, we describe some fundamental tools for dealing with pairs of random mutating strings. We first need to give a careful probabilistic exploration of random k-mer anchors on S = x₁x₂… and S′ = y_p+1y_p+2… . This requires a bit of work owing to the dependence between the random strings S, S′. For the rest of the paper, missing proofs can be found in the Supplemental Materials. Supplemental Figures S2 through S6 are visual aids for our proofs and can also be found in the Supplemental Materials.

We would like for M(i, j), M(i + 1, j + 1) to be independent, so finding the probability of k-mer matches is easy. However, in our model, it is not actually obvious a priori that M(i, j), M(i + 1, j + 1) are independent when j ≠ i. Consider M(1, 2) and M(2, 3). The former is a function of x₁, y₂; the latter, x₂, y₃. However, x₂ and y₂ are dependent in our model, so more work is needed. A convenient graphical representation of independence for the M random variables is the match graph, which we define below.

A match graph is shown in Supplemental Figure S1. The main reason for defining the match graph is the following theorem, which allows us to graphically determine independence of the A variables.

We will denote the random variables ${\sum }_{i=p+1}^{p+m}A\left(i,i\right)=$ $\sum _{i}A\left(i,i\right)=N_H\left(n,m\right)=N_H$ and ${\sum }_{i\ne j}A\left(i,j\right)=N_S\left(n,m\right)=N_S$, respectively, as the number of homologous anchors and spurious anchors; we will drop the dependence on n, m to simplify notation. These are key random variables that we wish to bound later on. The theorem below follows directly from Corollary 1 by linearity of expectation.

Bounding sums of k-mer random variables

We proceed to bound the distribution of the random variables N_S and N_H, which are sums of dependent random variables. We first bound N_S by computing the second moments and using variance-based bounds. To do this, we need to examine the independence structure of the A(i, j) random variables.

Intuitively, the first condition states that two anchors do not overlap too much; for example, A(1, 1) and A(2, 2) are not independent when k = 3. The intuition behind the second condition can be illustrated by the following situation in which k = 1 and θ ≈ 0: consider the anchors A(1, 5), A(5, 1). Because it is likely that x₁ = y₁ and x₅ = y₅, if x₁ = y₅, then x₅ = y₁ with high probability, so x₁ = y₅ is not independent of x₅ = y₁.

We can also bound expected values of $N_H^2,N_S^2$, and N_SN_H, giving us variance estimates.

Now we can use the variance bound and Chebyshev's inequality to get the result below. Note the bound uses k = C log n; we will prefer this form for quantities directly used for proving the main result. The label F1 in the theorem refers to the particular event space for which the bound always holds. We will label each proposition that holds with high probability with the event space that we are operating in. We will continue this convention for the rest of the paper as this will be useful when computing our final bounds.

If C > 3, then for large n, N_S = 0 with high probability, and our analysis would be easy. However, we want C as small as possible. It turns out we can make C ≈ 2 for reasonable θ, significantly tightening our bounds.

For N_H, we can get a stronger exponential bound because of its independence structure. A(i, i)s, which we call homologous anchors, are only dependent in a small neighborhood around i of size k because k-mers on nonoverlapping substrings are independent. This is called k-dependence (not to be confused with k-independence) and is used by Blanca et al. (2022) to show N_H is asymptotically normal. Concentration bounds can also be translated in the k-dependent scenario (Janson 2004).

Proof of nonsketched main result

To prove the main result on seed-chain-extend without sketching, we will need to bound three quantities in expectation: the runtime of chaining, the recoverability of chaining, and the runtime of extension.

We first bound O(E[Nlog N]), the expected runtime of chaining (and also anchor sorting). Note that $\mathbb{E}\left[N\right]\le {\displaystyle \frac{1}{n^{C-2}}+m{n}^{-C\alpha }.}$ We would like E[Nlog N] = O(E[N] log E[N]) = O(mn^−Cαlog m) to hold. Unfortunately, Jensen's inequality only gives E[Nlog N] ≥ E[N] log E[N] because xlog x is convex. However, with a bit more work, we get the following:

Now we bound the expected recoverability of the chain. Given S and S′, let a homologous gap of size ℓ + k − 1 bases be an interval of ℓ consecutive k-mers for which no homologous anchors exist (i.e., the k-mers are mutated). In the context of a chain, a homologous gap will refer to a gap flanked by two homologous anchors. Technically, if ℓ consecutive k-mers are uncovered but are flanked by two homologous anchors, this gives ℓ − k + 1 uncovered bases. We will ignore these factors of k as we will show that they are asymptotically small. It turns out homologous gaps grow relatively slowly in n with high probability.

In a chain, gaps may also be flanked by one or two spurious anchors. We call these nonhomologous gaps. We first bound break lengths, which will imply good recoverability, and bound nonhomologous gaps later on.

The idea behind proving the above propositions is to work in a space of events $\mathrm{F}1\cap \mathrm{F}2$_, where “bad events” do not occur and any optimal chain has good recoverability. Because this space of bad events is small, they do not contribute to our expected value too much. This finishes the recoverability result for the main theorem.

The last step is to bound extension running time, and this comes down to separately bounding the size of the homologous and nonhomologous gaps in any optimal chain. We bound the runtime of extension through homologous gaps by directly calculating the expectation through all possible homologous gaps. We then show that the runtime through nonhomologous gaps is small and does not contribute to the asymptotic term.

Sketching and local k-mer selection

Now consider not selecting all of the k-mers in a string, but only a subset of them during the initial seeding step. This allows one to chain only a subset of the k-mers, potentially providing runtime savings.

We use the open syncmer method (Edgar 2021). Given a string, we take all k-mers of the string and break the k-mer into s-mers with s < k. There are k − s + 1 s-mers in the k-mer. We select or seed the k-mer if the smallest s-mer (subject to some ordering, which we choose as uniform random) is in the $\lceil {\displaystyle \frac{k-s+1}{2}}\rceil$-th one-indexed position; we call a selected k-mer an open syncmer. Given repeated smallest s-mers, we take the rightmost one to be the smallest. Finding the smallest s-mer among the k − s + 1 s-mers in a k-mer takes k − s + 1 iterations, so finding all open syncmer seeds in S′ takes O((k − s + 1)m) = O(mk) = O(mlog n) time.

The expected fraction of selected k-mers over a string with i.i.d. uniform letters is called the density, and it is ${\displaystyle \frac{1}{k-s+1}}$ for the open syncmer method (up to a small error term $O\left({\displaystyle \frac{{\left(k-s+1\right)}^2}{{\sigma }^s}}\right)$ which we will ignore) (see (Zheng et al. 2020). We will let c be the reciprocal of the density, so c = (k − s + 1).

The original open syncmer definition had a parameter t, where a k-mer was selected if the smallest s-mer was in the tth position; we proved previously (Shaw and Yu 2022) that the optimal t is $\lceil {\displaystyle \frac{k-s+1}{2}}\rceil$ with respect to maximizing the number of conserved bases from k-mer matching. The reason we choose open syncmers is primarily to the following fact, which was shown by Edgar (2021):

Theorem 7 follows by examining the smallest s-mer in a k-mer and noticing that in the next overlapping k-mer, the locations for the new smallest s-mer are restricted. This theorem is the reason we use open syncmers and is crucial to our proofs. The spacing property makes selected open syncmers a polar set (Zheng et al. 2021); other methods also give rise to polar sets (Frith et al. 2021, 2023) but open syncmers seem to perform well empirically (Dutta et al. 2022; Shaw and Yu 2022; Frith et al. 2023) and are easy to describe. For the rest of the section, we will assume c = k − s + 1 is odd, so $\lceil {\displaystyle \frac{c}{2}}\rceil ={\displaystyle \frac{c+1}{2}}$.

Let A(i, j) be the random variables as defined before. Let J(i) be the indicator random variable set to 1 when the ith k-mer is selected on S as an open syncmer; J′(j), similarly for the jth k-mer on S′. We now wish to calculate E[J(i)J′(j)A(i, j)], the probability that a k-mer match exists and the k-mer is an open syncmer.

The above follows directly from Lemma 10. As expected, subsampling to a ${\displaystyle \frac{1}{c}}$ fraction of the k-mers gives ${\displaystyle \frac{1}{c}}$ expected hits. Note the important property of context independence used in the proofs: If A(i, j) = 1, then J(i) = J′(j). This property is not satisfied if one were sampling k-mers randomly or using minimizers (Roberts et al. 2004).

We now deduce the sketched moment bounds on N*:

The sketched moment bounds allow us to start reanalyzing the crucial propositions in the subsection “Proof of nonsketched main result” but in the context of sketching. The first main result is that sketching reduces the chaining time as one would expect.

The second result we wish to highlight is the new bound on extension runtime through homologous gaps.

It will turn out that in the sketched case, extension over homologous gaps also dominates runtime. Because we sketch with density ${\displaystyle \frac{1}{c}}$, this result states that we can sketch with decreasing density and have the same asymptotic extension time as without sketching, which is surprising. The crucial fact we use in the proof of Lemma 13 is Theorem 7, which shows that open syncmer seeding weakens the k-dependence of the seeds because they are now at least (c + 1)/2 bases apart. This tightens the bound in Theorem 7 and also shows that the sketched maximum homologous gap size is Θ (g(n)), with g(n) as in Lemma 6.

If we used other methods such as closed syncmers (Edgar 2021) or FracMinHash (Irber et al. 2022; Rahman Hera et al. 2023), then we would recover the O(c · mn^Cαlog n) result in Lemma 13 but not the $O\left({\displaystyle \frac{1}{c^2}m{n}^{C\alpha }{\log }^3\left(n\right)}\right)$ result. Thus, we have saved a log factor by using open syncmers when we let c = Θ (log n). The bulk of Section E in the Supplemental Materials is dedicated to proving Lemma 13.

The key quantities used for the rest of the proofs are $N_S^{\ast }$ and g′(n), the new maximum homologous gap size. We will prove g′(n) = Θ (g(n)) with g(n) as in Lemma 6, so the only asymptotic difference in these bounds is a factor of ${\displaystyle \frac{1}{\sqrt{c}}}$ in $N_S^{\ast }$. Because we want $N_S^{\ast }$ to be small, this does not affect downstream analysis. It follows that the proofs of the rest of the lemmas in the subsection “Proof of nonsketched main result” can be replicated almost verbatim. We give a sketch of these proofs in the Supplemental Materials.

Software availability

Scripts for regenerating the simulated experiments are available at GitHub (https://github.com/bluenote-1577/basic_seed_chainer/). The aligner used for the real nanopore experiments is available at GitHub (https://github.com/bluenote-1577/sce-aligner/) along with scripts for regenerating figures. All scripts are also available as Supplemenal Code.