Sampling k-mers is a ubiquitous task in sequence analysis algorithms. Sampling schemes such as the often-used random minimizer scheme are particularly appealing as they guarantee at least one k-mer is selected out of every w consecutive k-mers. Sampling fewer k-mers often leads to an increase in efficiency of downstream methods. Thus, developing schemes that have low density, i.e. have a small proportion of sampled k-mers, is an active area of research. After over a decade of consistent efforts in both decreasing the density of practical schemes and increasing the lower bound on the best possible density, there is still a large gap between the two.
We prove a near-tight lower bound on the density of forward sampling schemes, a class of schemes that generalizes minimizer schemes. For small w and k, we observe that our bound is tight when k≡1(modw). For large w and k, the bound can be approximated by 1w+k⌈w+kw⌉. Importantly, our lower bound implies that existing schemes are much closer to achieving optimal density than previously known. For example, with the current default minimap2 HiFi settings w=19 and k=19, we show that the best known scheme for these parameters, the double decycling-set-based minimizer of Pellow et al. is at most 3% denser than optimal, compared to the previous gap of at most 50%. Furthermore, when k≡1(modw) and the alphabet size σ goes to ∞, we show that mod-minimizers introduced by Groot Koerkamp and Pibiri achieve optimal density matching our lower bound.
Minimizer implementations: github.com/RagnarGrootKoerkamp/minimizers ILP and analysis: github.com/treangenlab/sampling-scheme-analysis.
For over a decade, k-mer sampling schemes have served as a ubiquitous first step in many classes of bioinformatics tasks. By sampling k-mers in a way which ensures that two similar sequences will have similar sets of sampled k-mers, sampling schemes enable methods to bypass the need to compare entire sequences at the base level and instead allow them to work more efficiently using the sampled k-mers.
Local sampling schemes satisfy a window guarantee that at least one k-mer is selected out of every window of w consecutive k-mers. Most schemes used in practice, such as the random minimizer scheme (Schleimer et al. 2003, Roberts et al. 2004), are forward schemes that additionally guarantee that k-mers are sampled in the order in which they appear in the original sequence. These properties are particularly appealing since they guarantee that no region is left unsampled.
As the purpose of these schemes is to reduce the computational burden of downstream methods while upholding the window guarantee, the primary goal of most new schemes is to minimize the density, i.e. the expected proportion of sampled k-mers. Over the past decade, many new schemes have been proposed that obtain significantly lower densities than the original random minimizer scheme.
For example, there are schemes based on hitting sets (Orenstein et al. 2016, Marc¸ais et al. 2017, 2018, DeBlasio et al. 2019, Ekim et al. 2020, Pellow et al. 2023, Golan et al. 2024), schemes that focus on sampling positions rather than k-mers (Loukides andPissis 2021, Loukides et al. 2023), schemes that use an ordering on t-mers (t < k) to decide which k-mer to sample (Zheng et al. 2020, Groot Koerkamp andPibiri 2024), and schemes that aim to minimize density on specific input sequences (Zheng et al. 2021b, Hoang et al. 2022). All of these improvements notwithstanding, it is still unknown how close these schemes are to achieving minimum density.
A trivial lower bound on density given by the window guarantee is 1 w , and recently Groot Koerkamp and Pibiri (2024) improved the bound of Marc¸ais et al. (2018) from
w þ k -0:5 . However, for many practical values of w and k, there is a sizeable gap between these lower bounds and the density of existing schemes. This raises the question whether schemes with density much closer to 1 w exist, but have not been found yet, or whether existing schemes are already very close to optimal and it is the lower bound that needs improvement. Our new lower bound closes most of the gap, and thus answers this question: Indeed, especially for k ≥ w, the best existing schemes have near-optimal density in many cases. This allows future research to focus on improving other sampling scheme metrics, such as the conservation described by Edgar (2021) and Shaw and Yu (2022).
Main lower bound theorem. We prove a novel lower bound on the density of forward schemes that is strictly tighter than all previously established lower bounds for all w, k, and alphabet size σ: Theorem 1. Let f be a (w, k)-forward sampling scheme and M σ ðpÞ count the number of aperiodic necklaces of length p over an alphabet of size σ. Then, the density of f is at least
where the middle inequality is strict for w > 1.
We prove that this bound can be extended to work for more general classes of sampling schemes, such as the local schemes described by Marc¸ais et al. (2018) and the multi-local schemes described by Kille et al. (2023).
Comparison with optimal schemes for small parameters. We show that our lower bound is tight for some small w, k, and σ by using an integer linear program to construct schemes whose density matches our lower bound. This marks the first time that there is an analytical description of a tight minimum density of any forward scheme. We conjecture that when k � 1ðmod wÞ, there exist schemes with density matching our lower bound.
Comparison with practical schemes for large parameters. To show that our bound is significantly closer to the density achieved by existing schemes compared to previous lower bounds, we replicate the benchmark from Groot Koerkamp and Pibiri (2024) for a selection of w and k (Fig. 3). For example, with the default minimap2 (Li 2018) HiFi settings w ¼ 19 and k ¼ 19, the lower bound goes up from 50% of the density achieved by the double decycling based method to 97% of the achieved density (Table 1).
Analysis of the mod-minimizer. Finally, our new lower bound implies that the mod-minimizer scheme (Groot Koerkamp and Pibiri 2024) is optimal when k � 1ðmod wÞ and σ is large. Indeed, for the ASCII alphabet (σ ¼ 256), the mod-minimizer scheme density is consistently within 1% of the lower bound when k � 1ðmod wÞ (Supplementary D, Fig. S4).
Notation. We begin by defining some necessary notation, as well as definitions of mathematical concepts that will be used throughout the work. We use ½n� to refer to the set f0; 1; . . . ; n -1g. The alphabet is denoted by Σ and has size σ :¼ jΣj, with σ ¼ 4 for DNA. The expression ajb indicates that a divides b. The summation P ajb is over all positive divisors a of b. We use a mod m for the remainder (in ½m�) of a after dividing by m and we use a � b ðmod mÞ to indicate that a and b have the same remainder modulo m. Given a string W, W½i; jÞ) refers to the substring of W containing the characters at 0-based positions i up to j -1 inclusive. For two strings X and Y, XY represents the concatenation of X and Y.
Classes of sampling schemes. There are multiple established classes of sampling schemes. We begin by drawing a distinction between schemes with and without a window guarantee that guarantees that at least one every w k-mers is sampled. While schemes without a window guarantee, such as fracminhash (Irber et al. 2022), are often efficient to compute, the lack of a guarantee on the distance between sampled k-mers makes them ineffective or inefficient for certain tasks such as indexing and alignment. Indeed, we only consider schemes with a window guarantee: Definition 1. A (w, k)-local scheme with window guarantee w and k-mer size k on an alphabet Σ corresponds to a sampling function f :
In other words, given a window of w þ k -1 characters (w consecutive k-mers), the output of the sampling function f (W) is an integer in ½w� which represents the index of the sampled k-mer in W. Recently, Kille et al. (2023) proposed a generalization of (w, k)-local schemes which samples at least s k-mers out of every w instead of at least 1 and we extend our results to these more general schemes in Supplementary A.
Local schemes have no restrictions on which of the w kmers can be selected for each window, but forward schemes are a subset of local schemes that enforce the restriction that they never select a k-mer which occurs before a previously selected k-mer. Koerkamp and Pibiri (2024) showed that 1:5 w þ k -0:5 is also a lower bound, 1 w is tighter for all of the parameter choices in the table. For SSHash (Pibiri 2022), we show parameters used for indexing a single human genome.
Definition 3. The density dðf Þ of a sampling scheme f is defined as the expected proportion of sampled positions from an infinite, uniformly random string.
For a further background on types of sampling schemes, we refer to Shaw and Yu (2022), Zheng et al. (2023), Groot Koerkamp and Pibiri (2024), and Ndiaye et al. (2024).
De Bruijn graphs. Let B n;σ ¼ ðV; EÞ denote the complete De Bruijn graph of order n, which has as vertices all strings of length n, V ¼ Σ n , and edges between vertices that overlap in n -1 positions, E ¼ fðX; X½1; nÞcÞjX 2 V; c 2 Σg. When σ is clear from the context or irrelevant for a particular discussion, it is omitted. It is worth noting that the vertices of
For each string s of length n, the n rotations of s induce a pure cycle in B n consisting of (up to) n vertices cyclically connected by edges. Note that when s is repetitive, e.g. a single repeated character or some other repeated string, the length of the cycle will be a divisor of n. These pure cycles are also called necklaces. The set of necklaces of length n corresponds to a partitioning of the vertices of B n into a vertex-disjoint set of pure cycles. We use C n to refer to this set of pure cycles of B n , and for c 2 C n , we write jcj for the number of vertices in the cycle.
When a string of length n has n unique rotations, the corresponding necklace is said to be aperiodic. The total number of necklaces and the number of aperiodic necklaces of length n are given by Moreau (1872) (and see also Riordan (1957)) as, respectively,
where φðpÞ is Euler's totient function that counts the number of integers in ½p� coprime to p. The formula M σ ðnÞ counting aperiodic necklaces follows from the formula for N σ ðnÞ via M€ obius inversion (M€ obius 1832), where μ is the M€ obius function defined to be 0 if n is divisible by a square (>1) and μðnÞ ¼ ð -1Þ q otherwise, where q is the number of prime factors of n. Charged contexts. The context of a window of length w þ k -1 in a sequence is the set of preceding windows that influences whether the current window samples a new position.
For a local scheme to select a new position, none of the previous w -1 windows may have selected the same k-mer as the current window. As a result, the context for local schemes consists of 2w þ k -2 characters: the current window of w kmers as well as the w -1 windows preceding the current window.
For a forward scheme, however, as soon as a window samples a different position than the preceding window, this position must be a new position. Thus, one needs only to consider the context of two consecutive windows of w kmers, for a total of w þ k characters.
When a sampling scheme selects a new position for the last window in a context, the context is charged. Marc¸ais et al. (2017) showed that the density of a scheme f can be defined as the proportion of contexts which are charged. In the case of forward schemes, each edge in B w þ k -1 represents a context, and the charged contexts are the edges (u, v) for which f ðuÞ 6 ¼ f ðvÞ þ 1.
Universal hitting sets. In 2021, Zheng et al. (2021a) related the density of forward and local schemes to the concept of universal hitting sets (UHS). A ðw; 'Þ-UHS is defined as a set of '-mers U such that any sequence of w adjacent '-mers must contain at least one '-mer from U. Theorem 1 of Zheng et al. (2021a) showed that when k ¼ 1, one can use the minimum size of a ðw; ' ¼ w þ kÞ-UHS to bound the density of ðw; k ¼ 1Þ-forward schemes, and the minimum size of a ðw; ' ¼ 2w þ k -2Þ-UHS to bound the density of a ðw; k ¼ 1Þ-local scheme.
In this section, we prove our main result: an improved lower bound on the density of forward sampling schemes. We first generalize some existing theorems to arbitrary w and k (Sections 3.1 and 3.2), after which our main theorem follows in Section 3.3.
We begin by considering a ðw ¼ 2; 'Þ-UHS. A ð2; 'Þ-UHS is equivalent to a vertex cover in B ' , i.e., a subset of vertices such that each edge in B ' is adjacent to at least one vertex in the subset. Lichiardopol (2006) used the fact that for every cycle C, at least djCj=2e of its vertices must be in a vertex cover, and obtained a lower bound on the size of a vertex cover by partitioning B ' into its pure cycles. We naturally extend this argument to obtain a lower bound on the cardinality of a ðw; 'Þ-UHS for any w ≥ 2. Notably, the pure cycle ð0011; 0110; 1100; 1001Þ has three vertices in the UHS, even though the lower bound given by Proposition 4 only requires it have 2. This is an example where the lower bound is not tight.
For certain values of w and ', such as when ' is prime or w ¼ 2 and ' is odd, Proposition 4 can be simplified to remove the summation and ceil function (Supplementary B).
Proposition 4 is the core of the proof of Theorem 1 and already has the right structure. The remainder of this section translates this result on universal hitting sets to a result on the density of sampling schemes.
Zheng et al. (2021a, Theorem 1) showed a connection between universal hitting sets and the density of sampling schemes when k ¼ 1. We naturally extend their result to k ≥ 1 for both local schemes (Lemma 5) and forward schemes (Corollary 6).
Lemma 5. Let f be a (w, k)-local scheme, and let C f be its corresponding set of charged contexts defined as the set of strings W of length 2w þ k -2 for which the last window W½w -
1ÞÞ not selected by any previous window:
Proof. For the sake of a contradiction, suppose there is a walk of length w in the De Bruijn graph of order ð2w þ k -2Þ, say ðW 0 ; . . . ; W w -1 Þ, that avoids cf Let S be the spelling of the walk, i.e., the sequence of length 3w þ k -3 such that S½i; i þ 2w þ k -2Þ ¼ W i . Since W w -1 6 2 C f and S contains W w -1 , this implies that on the last ðw þ k -1Þ-mer of W w -1 (i.e. S½2w -2; 3w þ k -3Þ), f selects an index j ≥ 2w -2 in S which has already been picked.
Since 0 ≤ f ð�Þ ≤ w -1 and j ≥ 2w -2, the first window that selects position j must begin at an index m ≥ w -1. Therefore, the context
as f selects a previously unselected position when applied to its last ðw þ k -1Þ-mer. By definition, W m -w þ 1 2 C f , contradicting our supposition and therefore C f is a ðw; 2w þ k -2Þ-UHS.
w Identically, one can consider contexts for a (w, k)-forward scheme f, which requires only verifying that the selection for a window of length w þ k -1 is distinct from the selection for the previous window. Therefore, the length of a context for forward f is only w þ k. As above, every w contexts must have at least one charged context, leading to the following conclusion: Corollary 6. If f is a (w, k)-forward scheme and C f is its corresponding set of charged contexts, defined as
As all contexts of a particular length ' are equally likely to occur in an infinite, uniform random string, the proportion of charged contexts corresponds to the density of the sampling scheme (Marc¸ais et al. 2017), i.e. dðf Þ ¼ jC f j=σ ' , where ' ¼ w þ k for forward schemes and ' ¼ 2w þ k -2 for local schemes. An example of the charged contexts of a (2, 2)-forward scheme and the corresponding UHS is depicted in Fig. 1.
We are now ready to state and prove our main theorem.
Theorem 1. Let f be a (w, k)-forward sampling scheme and M σ ðpÞ count the number of aperiodic necklaces of length p over an alphabet of size σ. Then, the density of f is at least
where the middle inequality is strict for w > 1.
Proof. Due to Corollary 6 and Marc¸ais et al. (2017), we can see that a (w, k)-forward sampling scheme of density dðf Þ implies a ðw;
By Proposition 4, this implies that every forward sampling scheme has a density of at least g σ ðw; kÞ, and hence dðf Þ ≥ g σ ðw; kÞ follows. For any p that divides w þ k, we haved p w e ≥ p w þ k d w þ k w e, with strict inequality when p ¼ 1 and w > 1. Substituting this in g σ ðw; kÞ, the middle inequality follows directly using the identity
that counts the number of strings of length w þ k partitioned by their shortest period. The last inequality follows directly from
As shown in Section 4, g σ ðw; kÞ is a tight bound for many small cases. Since its formula is somewhat unwieldy, 1 w þ k d w þ k w e can be used as an approximation that quickly approaches g σ ðw; kÞ (Fig. 2). Simple arithmetic shows that both g σ ðw; kÞ and 1 w þ k d w þ k w e improve the previous lower bound of 1:5 w þ k -0:5 of Groot Koerkamp and Pibiri (2024).
Given a (w, k)-local scheme f k , we can construct a ðw; k 0 ≥ kÞ-local scheme f k 0 of the same density by ignoring the last k 0 -k characters in each window, i.e. (Zheng et al. 2021a). It follows that the minimum density of a (w, k)-local or forward scheme is monotonically decreasing as k increases. However, as can be seen in Fig. 2, g σ ðw; kÞ is not a monotonically decreasing function. The local maxima appear to be at k � 1 ðmod wÞ, which motivates the following improved lower bound.
Theorem 2. For any (w, k)-forward scheme f, an improved lower bound g 0 is given by dðf Þ ≥ g 0 σ ðw; kÞ :¼ maxðg σ ðw; kÞ; g σ ðw; k 0 ÞÞ
where k 0 is the smallest integer ≥ k such that k 0 � 1 ðmod wÞ.
Remark 7. Similar to Theorem 1, Lemma 5 implies that any (w, k)-local scheme f has density at least dðf Þ ≥ g σ ðw; w þ k -2Þ. As this bound is in terms of g σ , the improved bound in Theorem 2 can be applied to local schemes as well, i.e., for any (w, k)-local scheme f, an improved lower bound is given by
4 Empirical tightness of our bounds
Here, we compare our bounds g σ and g 0 σ to existing lower bounds. Further, we show how tight these bounds are for small w, k, and σ by searching for optimal schemes via an integer linear programming (ILP) formulation. We also show how close existing sampling scheme densities are to g 0 σ for practical choices of w, k, and σ. Finally, we show when the recently described mod-minimizer scheme (Groot Koerkamp and Pibiri 2024) achieves optimal density as σ ! 1.
ILP description. We use an ILP to search for minimum density forward sampling schemes. In short, we use a single integer variable x W 2 ½w� for every window W of length w þ k -1 (corresponding to a vertex in B w þ k -1 ) that indicates the position of the chosen k-mer, and a single boolean variable y ðW;W 0 Þ for each edge in B w þ k -1 that indicates whether the corresponding context is charged. On each edge, we require that the scheme be forward. The objective function is to minimize the number of charged edges. To reduce the search space, we add an additional constraint corresponding to our lower bound g σ by requiring that for each pure Figure 2. Comparison of forward scheme lower bounds and optimal densities for small w, k, and σ. Optimal densities were obtained via the ILP and are plotted as black circles that are solid when the optimal density matches our lower bound, g 0 σ , and hollow otherwise. Each column corresponds to a parameter being fixed to the lowest non-trivial value, i.e., σ ¼ 2 in the first column, w ¼ 2 in the second column, and k ¼ 1 in the third column. Note that the x-axis in the third column corresponds to w, not k.
cycle of length jcj in B w þ k , at least djcj=we of the corresponding edges in B w þ k -1 are charged. Further details, including the ILP formulation for local schemes, can be found in Supplementary C.
Comparison against optimal schemes for small k. We used Gurobi (Gurobi Optimization, LLC 2024) to solve the ILP for all combinations of w, k, and σ such that 1 ≤ w ≤ 12; 1 ≤ k ≤ 12, and 2 ≤ σ ≤ 4 for both forward and local schemes and limited the runtime for each instance to 12 h on 128 threads. All results are reported in in Supplementary D, Table S2. While the additional constraint on pure cycles corresponding to g σ significantly sped up the search, for most large w, k, and σ, the ILP failed to terminate with an optimal solution in the allotted time. As a result, we restrict most of our analysis to the following three cases: fixed alphabet size σ ¼ 2, fixed window size w ¼ 2, and fixed k-mer size k ¼ 1 (Fig. 2).
For all ðw; k; σÞ where k � 1 ðmod wÞ (including when k ¼ 1), the minimum density exactly matches our lower bound g σ ðw; kÞ. Additionally, when σ ¼ 2 and w ¼ 2, the minimum density was equal to g 0 σ ðw; kÞ. Comparison against existing schemes for large k. Using a sequence of 10 million random characters over alphabet size σ ¼ 4, we approximated the density of recent sampling schemes using the benchmarking implementation from Groot Koerkamp and Pibiri (2024). To compare each density to the particular proportion of selected k-mers on a genomic sequence, we also ran all sampling schemes on the human Y chromosome (Rhie et al. 2023) after removing all non-ACTG characters. The densities of the best performing methods, Miniception (Zheng et al. 2020), double decycling-set-based minimizers (Pellow et al. 2023), and mod-minimizers (Groot Koerkamp and Pibiri 2024) are plotted in Fig. 3 along with random minimizers and lower bounds.
The ratio between the minimum achieved densities and lower bounds for a selection of (w, k) pairs used by existing k-mer-based methods are presented in Table 1. Additional results for σ 2 f2; 256g and w 2 f2; 50g are provided in in Supplementary D, Fig. S4.
The mod-minimizer has optimal density for large σ when w � k ðmod 1Þ. When w and k are constant and σ ! 1, the probability of duplicate characters in a window goes to 0. This implies that we can use t ¼ 1 for the mod-minimizer. When k � t ¼ 1 ðmod wÞ, the density of the mod-minimizer (Theorem 10 of Groot Koerkamp and Pibiri, 2024) is given by
The oð1='Þ term only accounts for duplicate t-mers, and hence disappears when σ ! 1. We get
Thus, the mod-minimizer has density equal to the lower bound provided by Theorem 1 when σ goes to 1 and w and k � 1 ðmod wÞ are fixed.
In practice, for σ ¼ 256 the mod-minimizer scheme is within 1% from optimal when k � 1 ðmod wÞ (Supplementary D, Fig. S4). When σ ¼ 4 (Fig. 3), a t > 1 must be used, causing the density plot to 'shift right' compared to the lower bound. Because of that, the modminimizer does not quite match the lower bound for practical values of σ.
Analytically, it is clear that g 0 σ ðw; kÞ is much larger than 1 w . In all cases, g 0 σ ðw; kÞ is nearly tight, if not completely. In particular, our bound is tight for all 40 tested parameter sets where k � 1 ðmod wÞ, leading us to our conjecture: Conjecture 1. For any w and k satisfying k � 1 ðmod wÞ, there exists a (w, k)-forward sampling scheme f such that dðf Þ ¼ g σ ðw; kÞ.
While the minimum size of a decycling set, i.e., a ðw ¼ 1; 'Þ-UHS, is well known to be N σ ð'Þ (Mykkeltveit 1972), very little is known about the minimum size of a ðw; 'Þ-UHS for finite w. In addition to providing the minimum density of a (w, k)-forward scheme for k � 1 ðmod wÞ, proving Conjecture 1 would also determine the minimum size of a ðw; ' ¼ w þ kÞ-UHS when k � 1 ðmod wÞ.
A natural investigation which follows our proposed lower bound is to determine the gap between g 0 σ ðw; kÞ and current forward scheme densities. Previously, the gap between known densities and lower bound was rather large, making it unclear how much more the density could be reduced.
In Table 1, we observe that existing schemes are already within 11% from the optimal density for practical values of w and k across different applications, and in many cases are even within 3% of the optimal density. In Fig. 3, we see that this difference holds not just for the specific (w, k) in Table 1, but for most k ≥ w. This is much more informative than the previous lower bound of 1=w, which implied that most current schemes are at most 50% denser than optimal for many of the parameters in Fig. 3.
For alphabets much larger than DNA ðσ ¼ 4Þ, such as the ASCII alphabet (σ ¼ 256), we observe that when k � 1 ðmod wÞ, the mod-minimizer scheme recently proposed by Groot Koerkamp and Pibiri (2024) is at most 1% denser than optimal and furthermore, we show that it is optimal as σ ! 1. This makes the mod-minimizer scheme the first practical scheme for which there exist finite parameters k and w for which it is close to optimal.
Our new bound for forward schemes always improves over 1=w and appears tight when k � 1 ðmod wÞ. This leads to an increasingly close bound for k 6 � 1 ðmod wÞ as k increases, but leaves a large gap when 1 < k < w. A better understanding of these small cases will be necessary to obtain a tight lower bound for all w and k. Based on Supplementary D, Figs. S3 and S4, one might conjecture that the double decyling-setbased methods of Pellow et al. (2023) are near-optimal, but subsequent work such as the greedy minimizer (Golan et al. 2024) has shown better schemes are possible. From Fig. 2, we already know that our lower bound is not always tight, so this leaves the question:
Open problem 1. How close can practical sampling schemes get to the density given by our lower bound?
For local schemes, though, our bound appears much less tight. We identified eight sets of ðw; k; σÞ where local schemes can obtain lower densities than their forward counterparts. In all cases, however, the difference between the local and forward densities was minuscule, with the largest difference of being found for ðw ¼ 4; k ¼ 2; σ ¼ 2Þ where the density decreased from 0.375 to 0.371 (Supplementary D, Table S2). Nevertheless, for some parameters, local schemes are able to achieve densities lower than our g 0 σ ðw; kÞ lower bound for forward schemes. Given the trend observed in Supplementary D, Table S2, we arrive at our final open problem:
Open problem 2. How much lower can the density of a (w, k)-local scheme be compared to a (w, k)-forward scheme?
Downloaded from https://academic.oup.com/bioinformatics/article/41/1/btae736/7922553 by guest on 01 May 2025
Kille et al.