{"paper_id":"015eef58-6e0a-4b4d-b71d-87b142bf7b23","body_text":"1\nFinimizers: Variable-length bounded-frequency\nminimizers for k-mer sets\nJarno N. Alanko, Elena Biagi, and Simon J. Puglisi\nAbstract—The minimizer of a k-mer is the smallest m-mer\ninside the k-mer according to some order relation < of the\nm-mers. Minimizers are often used as keys in hash tables in\nindexing tasks in metagenomics and pangenomics. The main\nweakness of minimizer-based indexing is the possibility of very\nfrequently occurring minimzers, which can slow query times\ndown significantly. Popular minimizer alignment tools employ\nvarious and often wild heuristics as workarounds, typically by\nignoring frequent minimizers or blacklisting commonly occurring\npatterns, to the detriment of other metrics (e.g., alignment recall,\nspace usage, or code complexity).\nIn this paper, we introduce frequency-bounded minimizers,\nwhich we call finimizers, for indexing sets of k-mers. The idea\nis to use an order relation < for minimizer comparison that\ndepends on the frequency of the minimizers within the indexed\nk-mers. With finimizers, the length m of the m-mers is not\nfixed, but is allowed to vary depending on the context, so that\nthe length can increase to bring the frequency down below a\nuser-specified threshold t. Setting a maximum frequency solves\nthe issue of very frequent minimizers and gives us a worst-\ncase guarantee for the query time. We show how to implement\na particular finimizer scheme efficiently using the Spectral\nBurrows-Wheeler Transform (SBWT) (Alanko et al., Proc. SIAM\nACDA, 2023) augmented with longest common suffix information.\nIn experiments, we explore in detail the special case in which\nwe set t = 1. This choice simplifies the index structure and\nmakes the scheme completely parameter-free apart from the\nchoice of k. A prototype implementation of this scheme exhibits\nk-mer localization times close to, and often faster than, state-\nof-the-art minimizer-based schemes. The code is available at\nhttps://github.com/ElenaBiagi/Finito.\nIndex Terms—Minimizer, Finimizer, Spectral Burrows-\nWheeler transform, pseudoalignment, k-mer, de Bruijn graph\nI. I NTRODUCTION\nI\nN recent years, use of minimizers [9], [35], [36], [40] has\nbecome widespread in bioinformatics, particularly as an\nalgorithmic tool to help manage large datasets for genome\nassembly [14] and for indexing tasks in pangenomics and\nmetagenomics [19], [25]–[27], [33]. The minimizer of a k-\nmer is the smallest m-mer inside the k-mer. The scheme is\nparameterized by the integers m ≤ k and an order relation <\nof the minimizers. Typically, k ≈ 30, m ≈ 18 and < is either\nthe lexicographic order, or an order determined by the hash\nvalues of the k-mers under some hash function.\nIn the context of indexing, a typical approach is to store\nlocation information for each minimizer of the data – for\nexample, the coordinates of their occurrences in the sequences\nbeing indexed, which may be, e.g., reads, assembled genomes,\nHelsinki Institute for Information Technology (HIIT) Department of Com-\nputer Science, University of Helsinki, Finland\nor unitigs 1 derived from such sequences. The index then\nconsists of a mapping of each minimizer to a list of their\noccurrences, along with the indexed sequences themselves.\nWith this information, it is possible to align query sequences\nof any length against the indexed data by computing the\nminimizers of the query sequence, locating them using the\nindex, and then checking whether the sequence around the\nminimizers matches (or approximately matches) the query.\nThe main weakness of this indexing strategy is that the\nquery time is proportional to the number of occurrences of\nthe query’s minimizers in the indexed data, which can be\nvery large. For example, in the unitigs of the human genome,\nminimizers can occur up to 3.6 ×104 times for k = 31\nand m = 20 [33]. Different tools attempt to get around this\nproblem in different ways. One possible method is to use\nmaximal exact matches as alignment seeds instead [24], but\nthis comes with its own drawbacks regarding performance and\nsensitivity. The most popular minimizer-based alignment tool,\nMINIMAP 2 [25], sidesteps the issue of frequent minimizers\nby ignoring minimizers that occur more often than a user-\ndefined threshold. The SSH ASH compact index for k-mer\nlookup [33] has a fallback for frequent minimizers, which can\nadd a significant amount of space to the basic index. Nyström-\nPersson et al. [30] observe that handling of minimizers in k-\nmer counting software can be particularly involved:\n“The k-mer counter KMC2 [12] introduced minimizer\nsignatures, which order m-mers lexicographically, except\nthat in order to reduce data skew, m-mers starting with\nAAA or ACA are given lower priority, and m-mers con-\ntaining AA anywhere are also avoided, except for AA at\nthe start . . . [t]his ordering is also used by Gerbil [15] and\nby FastKmer [17] in a modified form with some additional\nrules2”.\nAll these workarounds can have a detrimental effect on other\nimportant metrics, such as alignment recall, index size, and\ncode complexity. Another general difficulty with the use of\nminimizers is that the minimizer length m needs to be fixed\nbefore indexing, even though it may be hard for the user to\nknow what an efficacious value of m would be in advance.\nIn this paper, we introduce frequency-bounded minimizers,\nwhich we call finimizers. The idea is to use an order relation\n< for minimizer comparison that depends on the frequency\nof the minimizers themselves. With finimizers, the length m\nof the m-mers is not fixed, but is allowed to vary, so that the\n1Non-branching paths in the de Bruijn graph of the sequences.\n2Our emphasis.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n2\nlength can increase to bring the frequency down below a user-\ndefined maximum threshold t. Setting a maximum frequency\nthis way solves the issue of very frequent minimizers and gives\nus a worst-case guarantee for the query time.\nOur specific context throughout is the indexing of unitigs\n(or, more generally, any disjoint spectrum preserving string\nset, DSPSS). In this context, a set of finimizers is guaranteed\nto exist: when the input is a DSPSS, all k-mers occur only\nonce, so there are always low-frequency m-mers.\nThe principal advantages of finimizers for indexing are their\nbounded worst-case localization time and that they elide the\nneed for the user to fix the minimizer length m. While it\nmay appear we are exchanging one parameter for another (i.e.,\nminimizer length for finimizer frequency), we show that it is\npossible to efficiently find sets of finimizers with frequency\nt = 1, making some specific manifestations of finimizers\nessentially parameter free.\nIndeed, the main contribution of this article is to de-\nfine and explore such a finimizer scheme—Shortest Unique\nFinimizers—that is fast to compute, compact to store, and can\noutperform state-of-the-art minimizer-based schemes when\napplied to k-mer localization. In practice, our prototype imple-\nmentation of shortest-unique finimizers exhibits a k-mer local-\nization index that takes 25-30 bits per k-mer on a wide range\nof datasets, making it only marginally larger than SSH ASH.\nSimultaneously, it allows localization queries at comparable\nspeeds with asymptotic guarantees on performance and with-\nout the need to set parameters at construction time.\nAt face value, finimizers may seem to come with their own\nset of drawbacks. Firstly, unlike with minimizers, computation\nis no longer local, as an efficient means of determining the\nfrequency of a (potential) finimizer in the whole dataset is re-\nquired. Moreover, the fixed-length nature of minimizers means\nthey are particularly amenable to localization via sliding-\nwindow hashing schemes: the potential minimizer at a given\nposition is always the substring of length m starting there.\nWith their variable length, finimizers afford no such luxury.\nHowever, we show how both these problems can be elegantly\nand scalably solved via the use of the Spectral Burrows-\nWheeler transform (SBWT) [4] in lieu of hashing. As a bonus,\nwe also describe a fast heuristic for unitig orientation that\ncan significantly reduce the size of SBWT-based indexes —\nindeed, by 25% or more in our experiments.\nApart from efficient k-mer localization queries, we envisage\na vast array of possible future applications for finimizers. We\nlist some of these ideas at the end of the paper.\nRelated Work. Finimizers are tangentially related to,\nthough decided different from, frequency-aware minimizer\nschemes [9], [20]. Chikhi et al. [9] suggest ordering minimizer\ncandidates by frequency. In the robust winnowing scheme used\nby winnowmap [20] the set of minimizers occurring more\nthan 1024 times is stored and those m-mers are assigned a\nlower weight during minimizer selection. Another superficially\nrelated technique is universal hitting sets, where for integers\nk and L > k , a set of k-mers is a universal hitting set\nif every possible L-length sequence must contain a k-mer\nfrom the set [32]. We emphasise that, unlike finimizers,\nthese approaches still use fixed-length m-mers and provide\nno guarantees on minimizer frequencies.\nIn the context of k-mer localization queries, the state-of-the-\nart in minimizer-based schemes is the recent and fast S SHASH\nof Pibiri [33]. To give SSHASH its due, it does offer asymptotic\nguarantees on localization time; however it does so by making\na special case of frequent minimizers, an approach that does\nnot always pan out well. For example, on a set of Nanopore\nE.coli reads 3, setting minimizer length m = 14 increases\nrelative index size by almost 30% (from 12.7 to 16.2 bits/kmer)\neven though k decreases from 31 to 21.\nRoadmap. The next section lays down notation and basic con-\ncepts used throughout. Section III formally defines finimizers,\nincluding the shortest unique finimizer scheme. Section IV\nthen describes how to compute and query shortest unique fini-\nmizers as well as other details of our k-mer localization index.\nIn Section V, we describe a heuristic for orienting unitigs that\ncan significantly reduce SBWT-index size. In Section VI, we\nreport statistics on finimizer length and density, and the results\nof a localization query benchmark. Conclusions, reflections\nand potential future directions are then offered in Section VII.\nII. P RELIMINARIES\nThroughout this paper, a string X[1..n] is a sequence of |X|\nsymbols over the DNA alphabet Σ = {A, C, G, T}. The empty\nstring is denoted ϵ and |ϵ| = 0. The substring of X starting\nat symbol i and ending at symbol j is denoted X[i..j]. We\nalso use the half-open interval notations X(i..j] = X[i + 1..j]\nand X[i..j) = X[i..j −1]. A prefix is a substring with starting\nposition 1 and a suffix is a substring with ending position n.\nA k-mer refers to a (sub)string of length k. A minimizer is\nthe smallest m-mer of k-mer:\nDefinition 1. (Minimizer) Let Y be a k-mer, and m be an\ninteger such that m ≤ k. The minimizer of Y is the smallest\nm-mer of Y according to a given order relation of minimizers.\nThis order relation could be the simple lexicographic order or\nsomething more complex. The set of k-mers of a string X is\ncalled the k-spectrum of X.\nDefinition 2. (k-Spectrum) The k-spectrum Sk(X) of string\na X is the set of distinct k-mers {X[i..i + k − 1] | i =\n1, . . . ,|X| − k + 1} The k-spectrum Sk(X1, . . . , Xm) of a\nset of strings X1, . . . , Xm is the union Sm\ni=1 Sk(Xi)\nA padded k-spectrum adds a minimal set of $-padded dummy\nk-mers to ensure that in the de Bruijn graph (to be defined\nshortly), every non-dummy k-mer has an incoming path of\nlength at least k:\nDefinition 3. (Padded k-Spectrum) Let R = Sk(X1, . . . , Xm)\nbe a k-spectrum with alphabet Σ, and let R′ ⊆ S be the set\nof k-mers Y such that Y [1..k − 1] is not a suffix of any k-mer\nin R. The padded k-spectrum is the set S+\nk (X1, . . . , Xm) =\nR∪{$k}∪S\nY ∈R′ {$k−iY [1..i] | i = 1, . . . , k−1}, with special\ncharacter $ ̸∈ Σ and $ < c for all c ∈ Σ.\n3SRA identifier SRR25689478, see also Section VI-B.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n3\nFor example, if X1 = ACGT, X2 = GACG and k = 3,\nthen S3(X1, X2) = {ACG, CGT, GAC} and S+\n3 (X1, X2) =\n{ACG, CGT, GAC, $$$, $$G, $GA}. See Figure 3 for a more\nelaborate example.\nThe de Bruijn graph is an important concept in k-mer in-\ndexing. In this paper, we use a padded node-centric definition,\nwhere the de Bruijn graph of a set of strings X1, . . . , Xm is an\nedge-labeled graph G = (V, E), where V = S+\nk (X1, . . . , Xm)\nand there is an edge from k-mer u to v iff u(1..k] = v[1..k).\nSee Figure 1 for an example. The label of an edge e = (u, v),\ndenoted with ℓ(e) or ℓ(u, v), is the character v[k] (with this,\n$ never appears as an edge label). The label of a path of\nedges is the concatenation of the edge labels in the order\nof the path, and the label of a path of nodes v1, . . . , vm is\nv1 · ℓ(v1, v2) · · · ℓ(vm−1, vm), with possible dollar symbols in\nv1 removed.\nA unitig is a maximal simple node path v1, . . . , vm in a\nde Bruijn graph such that outdegree(v i) = 1 for i > 1 and\nindegree(vi) = 1 if i < m . Maximality in this context means\nthat the path can not be extended in either direction without\nrepeating nodes or breaking indegree or outdegree constraints.\nThe set of all unitigs covers all nodes of the graph and is\nuniquely defined if a canonical choice of starting point is made\nfor disjoint cycles (for example, the lexicographic minimum).\nDefinition 4. (Disjoint Spectrum-Preserving String Set\n(DSPSS)) A set of strings X1, . . . , Xm is a disjoint spectrum-\npreserving string set of order k of a set of strings Y1, . . . Yn\niff Sk(X1, . . . , Xm) = Sk(Y1, . . . , Yn) and all k-mers in\nX1, . . . , Xm are distinct.\nThe set of node path labels of all unitigs in the de Bruijn\ngraph of a set of strings in an example of a DSPSS for that\nset of strings [34], and the smallest DSPSS can be computed\nin linear time using the Eulertigs algorithm [38].\nThe main operation we are interested in supporting on a\nk-spectrum is the localization query:\nDefinition 5. (k-mer localization) Given a disjoint spectrum-\npreserving string set D and k-mer Y , localize(Y ) returns\neither ⊥ if the k-mer does not occur in any X ∈ D , or\notherwise a pair (i, j), where i is the index of the string\nXi ∈ D that contains Y , and j is the starting point of Y\nin Xi.\nGiven a DSPSS D = X1, . . . , Xm, it will be important to\nmake a distinction between DSPSS occurrences and DBG\noccurrences of substrings:\nDefinition 6. (DSPSS occurrence). Given a DSPSS D =\nX1, . . . , Xm, a DSPSS occurrence of string Y is a pair (i, j)\nsuch that Xi[j..j + |Y | − 1] = Y .\nDefinition 7. (DSPSS frequency) Given a DSPSS D =\nX1, . . . , Xm, the DSPSS frequency of Y is the number of\ndistinct DSPSS-occurrences of Y .\nDefinition 8. (DBG occurrence) Given a DSPSS D =\nX1, . . . , Xm, a DBG occurrence of Y is a node of the de\nBruijn graph of D that has an incoming path with edges\nspelling Y .\nDefinition 9. (DBG frequency) Given a DSPSS D =\nX1, . . . , Xm, the DBG frequency of Y is the number of DBG\noccurrences of Y .\nNote that the number of times a substring occurs in the\nde Bruijn graph is not necessarily equal to the number of\noccurrences in the DSPSS, since, for example, the labels of\nneighboring unitigs in the unitig DSPSS overlap by k − 1\ncharacters. This is illustrated in Figure 2, where the DSPSS\nfrequency of CA is 5, but the DBG frequency is only 1. The\nDBG frequency equals the number of k-mers in the padded\nk-spectrum that have Y as a suffix.\nA path v1, . . . vn in a graph is considered branching if there\nexists i > 1 such that the indegree of vi is at least 2, or there\nexists i < n such that the outdegree of vi is at least 2. A\nstring Y is considered branching in a DBG if there exists a\nbranching path labeled with Y .\nA key tool in succinct data structures is support for\nrank and select queries on a bit string (or bit vector) Z\nof length n. For i ≤ n and z ∈ {0, 1}, rankz(i) is\nnumber of z’s in the first i bits of Z, and selectz(i) is the\nposition of the ith z in Z. Classical techniques need n+ o(n)\nbits to support each operation in O(1) time [29].\nIII. F INIMIZER FUNCTIONS\nDefinition 10. (Finimizer function) Let R be a set of k-mers\nand t be a positive integer. A function f : Σ k → Σ∗ is a t-\nfinimizer function for R if for all X ∈ Σk, f(X) is a substring\nof X, and f(X) has at most t occurrences in the de Bruijn\ngraph of R. The set {f(X) | X ∈ R} is called a t-finimizer\nset for R.\nA t-finimizer function always exists because the identity\nfunction fulfills the requirements: the number of occurrences\nof every full k-mer is 1, so the frequency requirements are\nfulfilled. For the purposes of indexing, however, we desire a\nfinimizer function such that the size of the image {f(X) | X ∈\nR} is small, and such that adjacent k-mers in the de Bruijn\ngraph likely have the same finimizer. Additionally, the function\nshould allow a compact representation that also ensures fast\nevaluation.\nThere are potentially many finimizer functions with these\nproperties. The one we focus on for the remainder of this\npaper is the function that picks the shortest substring with\nfrequency at most t.\nDefinition 11. (Shortest t-Finimizers) Let G be the de Bruijn\ngraph of a set of k-mers R, t ≥ 1 be an integer, X be a\nk-mer, and Y be a substring of X. We say Y is a shortest\nt-finimizer of X with respect to the k-mer set R if Y has at\nmost t occurrences in G and there does not exist a shorter\nsubstring of X with at most t occurrences in G.\nIf multiple shortest finimizer candidates exist, a consistent\ntiebreaker such as colexicographic comparison can be used\nto always pick the same finimizer.\nWhen t = 1, we call the finimizers Shortest Unique\nFinimizers (SUFs). As we shall see in Section IV, when\nimplemented with care, a finimizer function for shortest unique\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n4\nTA AAG TC CA GGT GA\n1 GTAA 3 AAGT 4 AGTC 15 ACAG 25 AGGT 9 AGGA\n2 TAAG 5 GTCT 16 CAGG 10 GGAA\n26 GGTA 11 GAAA\n20 GTAG\n21 TAGG\nCAGG\nAGGA\nAGGT\nGGAA\nAAGT\nGTAA\nGGTA\nGTAG\nTAGG\n$$$A\n$$$$\n$$AC\n$ACA\nACAG\nAGTC\nTAAG\nGTCTGAAA\nFig. 1. Left: A set of shortest unique finimizers of the set of unitigs U = {GT AAGTCT, AGGAAA, ACAGG,\nGT AGG, AGGT A}, k = 4. From the SFSk,1 of each unitig, we can obtain the set of shortest unique finimizers (SUFs)\nfor each k-mer. SFSk,1 = [⊥, ⊥, (2, 1, 8), (3, 1, 5), (3, 1, 11), (3, 1, 17), (2, 1, 10), (2, 1, 16)], [⊥, ⊥, ⊥, (2, 1, 7), (3, 1, 4), (3, 1, 3)],\n[⊥, (2, 1, 9), (2, 1, 6), (3, 1, 12), (4, 1, 14)], [⊥, ⊥, (2, 1, 8), (3, 1, 13), (4, 1, 15)], [⊥, ⊥, ⊥, (3, 1, 18), (2, 1, 8)]. These five strings, composed\nof 14 k-mers, can be represented with six finimizers. For the k-mer ACAG we choose CA over AC as CA<colexAC. The SUF for the k-mer\nGGT Ais T A, since it is shorter than GGT. Right: The de Bruijn graph of U. Red dashed edges are pruned from the graph because nodes\nthey point to can be reached from other (black) edges.\nAATACAG\nACAGT\nACAGA\nCAGTC\nCAGTA\nT\nA \nC\nA\nG A C A T A A \nFig. 2. A finimizer that is unique in the de Bruijn graph (right, k = 5) may not be unique in the unitig DSPSS (left). In our index structure,\nthe unitig occurrence that represents the de Bruijn graph occurrence is the one which is observed farthest from beginning of the unitig.\nfinimizers can be constructed in time linear in the total length\nof the input DSPSS, and evaluated for every k-mer of a query\nsequence Q in time O(|Q|).\nIV. C OMPUTING AND QUERYING SHORTEST UNIQUE\nFINIMIZERS\nIn this section we show how to index a set of unitigs for\nk-mer localization queries using Shortest Unique Finimizers\n(SUFs). The problem is related to the active line of research\non shortest unique substring queries [1]. Here, we present a\nsolution specialized to sets of k-mers by using the Spectral\nBurrows-Wheeler transform [4]. The methods generalize to\nany DSPSS, but here we focus on unitigs for concreteness.\nWe first show how the shortest unique finimizer for each k-\nmer in the DSPSS can be identified in linear time in the\ntotal length of the DSPSS using the Spectral Burrows-Wheeler\ntransform (SBWT) augmented with the longest-common-suffix\n(LCS) array, two data structures that we define below. See\nFigure 3 for an illustration of these structures.\nDefinition 12. (Spectral Burrows-Wheeler Transform\n(SBWT) [4]) Let R+ be a padded k-spectrum and let\nX1, . . . , X|R| be the elements of R+ in colexicographic\norder. The SBWT is the sequence of sets A1, . . . , A|R| with\nAi ⊆ Σ such that Ai = ∅ if i > 1 and Xi[2..k] = Xi−1[2..k],\notherwise Ai = {c ∈ Σ | Xi[2..k]c ∈ R+}.\nIt is helpful to think of the SBWT as encoding the set R\nof k-mers as a list arranged in colexicographical order. In the\ncontext of k-mer lookup, the SBWT allows navigation of that\nordered list via the operation ExtendRight, defined as follows.\nLet [s, e]α denote the colexicographic interval of string α,\nwhere s and e are respectively the colexicographic ranks of\nthe smallest and largest k-mer in the SBWT having substring\nα as a suffix. The right extension of the interval [s, e]α with a\ncharacter c ∈ Σ, denoted ExtendRight([s, e]α, c), is the interval\n[s′, e′]αc, or ⊥ if no such interval exists. ExtendRight can be\nanswered in O(1) time using a data structure of n|Σ|+o(n|Σ|)\nbits [4].\nIt is easy to see that an existential k-mer lookup algorithm,\nwhich tells if a query k-mer r exists (or not) in S, comes from\nthe repeated application of ExtendRight to the symbols of the\nquery k-mer, starting with the interval [1, n]ϵ, where ϵ is the\nempty string. Indeed the algorithm visits all intervals of the\ncolexicographically sorted k-mer spectrum that correspond to\nall prefixes of the query k-mer that are suffixes of k-mers in\nS. This also gives us the DBG frequency of every prefix of the\nquery k-mer since the DBG frequency of a string is equal to\nthe length of its colexicographic interval. To support efficient\nstreaming queries, we pair the SBWT with the LCS array:\nDefinition 13. (Longest Common Suffix (LCS) Array [2]) Let\nR+ be a colexicographically ordered padded k-spectrum and\nlet Xi be the colexicographically i-th k-mer of R+. The LCS\narray is an array of length |R+| such that LCS[1] = 0 , and\nfor i > 1, LCS[i] is the length of the longest common suffix of\nXi−1 and Xi.\nIn the above definition, we consider the empty string as\na common suffix of any two k-mers, such that the longest\ncommon suffix is well-defined for any pair of k-mers.\nThe LCS array can be constructed from the SBWT in O(n)\ntime [2]. It allows a different navigation operation—called\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n5\nContractLeft—on the colexicographically-ordered padded k-\nmer spectrum. In particular, operation ContractLeft([s, e]α) for\n|α| > 0 returns the interval [s′, e′]α[2..|α|], A ContractLeft\nquery can be answered by setting s′ to be the largest position\nsmaller or equal to s in LCS, such that LCS[s′] < |α| − 1 (a\nprevious-smaller-value query, PSV-query for short), or to 1 if\nno such position exists. Symmetrically, e′ can be found with\na next-smaller-value (NSV) query respect to e (see Figure 3).\nThere exist data structures that can answer both PSV and NSV\nqueries in constant time with modest space overhead [7], [18].\nIn practice however, we implement ContractLeft([s, e]α) by\nsimply scanning the LCS array left from s and right from e.\nThis works well because in our algorithms, in practice, the vast\nmajority of scans tend to be very short (at most 6 elements\non average in our experiments). The scanning solution is also\nmemory-local and is thus likely faster than more sophisticated\ndata structures with worst-case guarantees.\nThe ContractLeft operation in combination with the\nExtendRight operation allows us to compute in linear time\nthe k-bounded matching statistics:\nDefinition 14. (k-bounded matching statistics) The k-bounded\nmatching statistics of a query sequence Q[1..m], versus a\npadded k-spectrum R+, is the array MS k[1..m] such that\nMSk[i] = ( di, pi, ℓi), where di is the largest non-negative\ninteger such that Q(i − di..i] is a suffix of at least one k-\nmer in R+, pi is the DBG frequency of Q(i − di..i] and ℓi is\nthe start of the colexicographic interval of Q(i − di..i].\nWe can compute MSk efficiently by using a k-bounded variant\nof the BWT-based unbounded matching statistics algorithm of\nOhlebusch et al. [31]. The pseudocode is given in Algorithm\n1, and the proof of correctness is given in the Supplement. We\nwill use the MSk array during localization queries to find the\ncolexicographic ranks of all k-mers of the query present in\nthe index: all the k-mers ending at positions where the match\nlength is equal to k.\nBuilding on the idea of Algorithm 1, we can also compute\nthe shortest frequency-bounded suffixes:\nDefinition 15. (Shortest frequency-bounded suffixes) The\nshortest suffix statistics of a query sequence Q[1..m], versus\na padded k-spectrum R+, is the array SFS k,t[1..m] such that\nSFSk,t[i] = (di, pi, ℓi), where di is the smallest positive integer\nsuch that Q(i − di..i] is a suffix of pi (pi ≤ t) k-mers in R+,\nand ℓi is the start of the colexicographic interval ofQ(i−di..i].\nIf no such di exists, SFS k,t[i] = ⊥.\nThe SFSk,t array will be useful for identifying the shortest\nfrequency-bounded substring of every k-mer in a query. Com-\nputation of the SFSk,t array with unbounded k has been stud-\nied by Belazzougui and Cunial [5]. Their algorithm uses data\nstructures on the suffix tree topology of the data to implement\nleft contractions. We avoid most of the complications in their\nalgorithm by replacing the suffix tree topology with the LCS\narray, which uses only ⌊log(k−1)⌋+1 bits per k-mer since the\nLCS values are in the range [0, k − 1]. Algorithm 2 gives the\npseudocode in terms of left contractions and right extensions.\nWe give a sketch of the proof of correctness here. The full\ndetails can be found in the supplement.\nAlgorithm 1 k-bounded matching statistics.\nInput: A query string Q and a data structure supporting ExtendRight\nand ContractLeft on the padded k-mer spectrum.\nOutput: MSk[1..|Q|]\n1: [ℓ, r] ← [1, n] ▷ The current colexicographic interval\n2: d ← 0 ▷ The length of the current match, up to k\n3: MSk[1..|Q|] ← Array of length |Q|.\n4: for i = 1..|Q| do\n5: while d > 0 and ExtendRight([ℓ, r], Q[i]) = ∅ do\n6: [ℓ, r] ← ContractLeft([ℓ, r])\n7: d ← d − 1\n8: if ExtendRight([ℓ, r], Q[i]) ̸= ∅ then\n9: [ℓ, r] ← ExtendRight([ℓ, r], Q[i])\n10: d ← min(k, d + 1)\n11: MSk[i] ← (d, r − ℓ + 1, ℓ) ▷ Len, Freq, Colex start\nAlgorithm 2 Shortest frequency-bounded suffix statistics.\nInput: A query string Q, frequency upper bound t, and a data\nstructure supporting ExtendRight and ContractLeft on the padded\nk-mer spectrum.\nOutput: SFSk,t[1..|Q|]\nAssumes all characters of the query Q exist in the index at least\nonce.\n1: [ℓ, r] ← [1, n] ▷ The current colexicographic interval\n2: d ← 0 ▷ The length of the current match, up to k\n3: SFSk,t[1..|Q|] ← Array of length |Q| initialized to values ⊥\n4: for i = 1..|Q| do\n5: while ExtendRight([ℓ, r], Q[i]) = ∅ do\n6: [ℓ, r] ← ContractLeft([ℓ, r])\n7: d ← d − 1\n8: [ℓ, r] ← ExtendRight([ℓ, r], Q[i])\n9: d ← min(k, d + 1)\n10: while |ContractLeft([ℓ, r], 1)| ≤ t do\n11: [ℓ, r] ← ContractLeft([ℓ, r])\n12: d ← d − 1\n13: if r − ℓ + 1 ≤ t then\n14: SFSk,t[i] ← (d, r − ℓ + 1, ℓ)\nTheorem 1. Given constant-time ContractLeft and\nExtendRight queries on padded k-spectrum R+, Algorithm 2\ncomputes SFS k,t for query Q in time O(|Q|).\nProof. The algorithm maintains an invariant that at the end of\niteration i, the interval [ℓ, r] is the interval of the shortest suffix\nof Q[1..i] occurring as a suffix of at most t k-mers, if it exists,\nor otherwise [ℓ, r] is the longest suffix of Q[1..i] occurring in\nthe k-mers. Time is linear in |Q| since each query position is\nextended to and contracted from only once.\nA. Shortest Unique Finimizer Index Structure\nOur index structure is composed of four main components.\nWe describe these components in this subsection to give an\noverall view of our data structure. Subsequent sections show\nhow they are constructed efficiently and finally how they are\nused for querying.\nLet us denote the set of unitigs as U = {u1, u2, . . . , u|U |}\nsorted by the colexicographical order of their first k-mers. Let\nn denote the number of k-mers in the padded k-mer spectrum\nR+ of a set of U. It will be convenient for us to refer to\nk-mers in R+ by their colexicographical rank. Let Xi denote\nthe i-th k-mer in R+ in colexicographic order. Our SUF index\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n6\nfor U consists of the following components (see Figure 3 for\nan example):\n1. The k-mer index. The k-mer index includes the SBWT\nand the LCS array, which are used to identify the shortest\nunique finimizers using ExtendRight and ContractLeft opera-\ntions.\n2. The unitig index. The purpose of this component is to\nprovide random access to the unitig strings. It consists of the\nconcatenation U ′ = u1 ◦ . . . ◦ un, and the array E storing the\nendpoints of each unitig in U ′ in ascending order, compressed\nwith the Elias-Fano method (see [29]), with support for access\nand predecessor queries.\n3. The finimizer dictionary. This component is used to map\nfinimizers identified by their colexicographic ranks to their\nlocations (\"global offsets\") in U ′. The finimizer dictionary\nstores a bit vector F [1..n] such that F [i] = 1 iff Xi has a\nSUF as a suffix. Note that a bit vector suffices for F precisely\nbecause the SUFs are unique w.r.t the de Bruijn graph. The\nbit vector F is processed for constant-time rank queries. The\nfinimizer index also stores an array G such that G[i] is the\nending position in U ′ of an occurrence of the finimizer with\ncolexicographic rank i. Note that even though finimizers are\nunique in the DBG, they may occur in multiple unitigs strings.\nTo resolve this ambiguity, we always pick the uniquely defined\nunitig occurrence that occurs furthest from the start of a unitig.\nThe correctness of the query algorithm will depend crucially\non this choice.\n4. The branch dictionary. The purpose of this component\nis to identify k-mers that are the first k-mers of their unitigs.\nTo this end, we store a bit vector UStart[1..n] that marks the\nk-mers that are prefixes of unitigs. That is, UStart[i] = 1 iff\nthe Xi is a prefix of some uj ∈ U. The bit vector UStart is\npreprocessed for constant-time rank and select queries. This\ncomponent will be important when localizing k-mers that cross\nbranches in the DBG.\nB. Shortest Unique Finimizer Index Construction\nIn this section, we describe the construction of our index\nfor shortest unique finimizers. We begin by sorting the unitigs\nby the colexicographic ranks of their first k-mers, and con-\ncatenating the sorted unitigs to construct U ′ and E. Then,\nwe construct the SBWT using the algorithm described in [4]\nand from there, build the LCS array using the techniques in\n[2]. The bit vector U start is easily constructed by searching\nthe first k-mer of each unitig on the SBWT and marking the\ncolexicographic positions of those.\nAlgorithm 3 gives pseudocode for constructing the fin-\nimizer dictionary, consisting of the global offsets array G\nand the finimizer colexicographic marks bit vector F . The\nalgorithm processes the unitig strings one by one in the same\norder as they are concatenated in U ′. For each unitig, we run\nAlgorithm 2 to compute the shortest unique suffixes array,\nSFSk,1, (Def. 15) of the unitig. Then, we slide a window of\nlength k over the SFSk,1 array while maintaining a monotonic\nAlgorithm 3 Shortest Unique Finimizers dictionary construction\nInput: Sorted unitigs U = {u1, u2, . . . , u|U |}, Endpoints array\nE[1..|U |], SBWT(U ), LCS array of padded the k-spectrum of U.\nOutput: finimizer dictionary (G, F )\nG ← Empty list\nF, D, A ← Arrays of length n initialized to zeros\nMQ ← Monotonic queue of triples (len, ℓ, end)\nfor i = 1 . . . |U | do\nSUSk ← Run Algorithm 2 on ui with t = 1\nMQ.push(SUSk[1 . . . k − 1])\nfor kend = k . . . |ui| do\nMQ.push(SUSk[kend])\n(lenv, ℓv, endv) ← MQ.min()\nwhile endv − lenv < k end − k do\nMQ.pop() ▷ Ensure SUF is in the current k-mer\n(lenv, ℓv, endv) ← MQ.min()\nif D[ℓv] ≤ endv then ▷ Store only the rightmost\nA[ℓv] ← endv + E[i − 1] ▷ Edge case E[0] = 0\nF [ℓv] ← 1\nD[ℓv] ← endv\nif |MQ| > k then\nMQ.pop()\nfor j = 1 . . . |F | do\nif F [j] = 1 then G.append(A[j])\nreturn (G, F )\nqueue data structure4 to be able to identify the shortest unique\nsuffix inside the current k-mer window, with colexicographic\ntiebreaking. This process gives us the colexicographic rank and\nending position of the finimizer of each k-mer of the unitig.\nWe then mark these colexicographic ranks in the bit vector F ,\nand store corresponding the global offsets. If a global offset is\nalready stored for some colexicographic position, we overwrite\nit iff the finimizer ending position is further from the start of\nthe unitig than the currently stored offset. The array D in the\npseudocode keeps track of the distance of the stored global\noffset from the start of a unitig. For example, in Figure 2\nthe final global offset stored for CA would be the one that is\ncontained in the leftmost unitig of the picture.\nSee Figure 1 for an example of Shortest Unique Finimizers\nand Shortest Unique Suffixes array of a set of unitigs, and\nFigure 3 for the full index for these finimizers.\nC. Query Algorithm\nWe are now ready to describe our algorithm for localizing\nall k-mers of a query string Q in the set of unitigs, using\nthe shortest unique finimizer index described in the previous\nsection.\nThe pseudocode is given in Algorithm 4. The first step of\nthe query algorithm involves running Algorithm 1 to obtain k-\nbounded matching statistics in order to identify all k-mers of\nthe query that exist in the index. This returns an array, MSk,\nof size |Q|, where each element is a triple (len, freq, ℓ ),\n4A monotonic queue is a data structure that can report the minimum element\nin a variable-length sliding window over a stream of data x1, x2, . . . xn. The\nqueue has three operations: Push(x): Appends element x to the window. Min():\nRetrieves the minimum element x in the current window along with its index\ni in the window. Pop(): Removes the oldest element in the window. In our\nuse case, the elements are triples of integers (length, frequency, colex rank),\nordered lexicographically.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n7\nAlgorithm 4 Streaming Localization Queries\nInput: A query string Q.\nOutput: For each k-mer of Q that is present in the unitigs, prints\nthe id of the unitig that contains it, and the position of the k-mer\nwithin that unitig.\nMSk ← Run Algorithm 1 on Q\nSUSk ← Run Algorithm 2 on Q with t = 1\nprev ← −∞ ▷ \"Previous branch\" end in Q\nfor i = 1 . . . |Q| do\n(len, freq, ℓ) ← MSk[i]\nif freq = 1 and UStart[ℓ] = 1 then\nprev ← i ▷ Update rightmost branch\nelse if freq = ⊥ then\nprev ← −∞\nif len = k then ▷ Full k-mer match ending at i\nif prev > 0 then ▷ Rightmost branch is still valid\nunitig_id ← UStart.rank(MSk[prev].ℓ)\nl_off ← i − prev + 1\nelse ▷ No branches crossed\n(lenv, ℓv, endv) ← finimizer of Q[i − k + 1 . . . i]\ng_off ← G[F.rank(ℓv)]\nunitig_id ← 1 + index of predecessor of g_off in E\nl_off ← g_off −E[unitig_id−1] + (i − endv) − k + 1\nprint (unitig_id, l_off )\nlen is the length of the longest match up to that position,\nfreq is the frequency of that substring and ℓ is the start of\nits colexicographic interval in the SBWT. If len = k, the\nwhole kmer at position i is considered to be found and we\nwill proceed with its localization. As a second step, we run\nAlgorithm 2 and use a sliding window minimum algorithm\nto pick the finimizer for each k-mer that we know exists in\nthe index. This procedure is identical to the one described\nfor identifying SUFs in the unitigs to build the finimizer\ndictionary.\nWe can locate each found k-mer Y in the input unitigs in\ntwo different ways, using the branch dictionary or the finimizer\ndictionary, depending on whether a branch has been observed\nbefore the k-mer endpoint. For each position i of Q with\nMSk[i].freq = 1 the algorithm checks if the identified unique\nmatch is the suffix of the first k-mer of a unitig, and if so it\nassigns the value i to prev. The starting k-mers are marked\nby the U start bit vector, thus prev is assigned to i if we\nobserved a 1 in UStart at the colexicographic index stored in\nMSk[i].\nAs a consequence of this, if a branch has been observed\nprev will have a value bigger than 0, the offset of the\nclosest previous branch. In this case, the unitig identifier\nof Y is obtained, using the branch dictionary, with a rank\nquery on U start since unitigs have been sorted based on the\ncolexicographic rank of their first k-mers. The local offset\nof the k-mer Y is given by shifting the current position i\naccording to the offset of the observed branch, prev.\nIn the case in which no branch has been observed, we\nuse the finimizer dictionary. Given the found SUF v, at\ncolexicographic rank ℓv, we retrieve, in the global offset array\nG, the global offset (variable g_off in the pseudocode) of the\nunique k-mer suffixed by the finimizer corresponding to ℓv,\nwhich which is at G[F.rank(ℓv)]. The unitig identifier of Y\ncan be found with a predecessor query on the unitigs endpoints\narray E.pred(g_off ). The local offset of the k-mer is calculated\nby subtracting from g_off the cumulative length of the previous\nunitig in the concatenation and shifting the resulting position\naccording to the offset of the finimizer in Y .\nExamples covering both of the main cases in the localization\nalgorithm can be found in Figure 3. A full proof of the\ncorrectness is included in the supplementary material.\nD. Implementation Details\nThis search algorithm can be improved with what we call a\nunitig walk after each found k-mer. This consists of a bitwise\ncomparison of the characters following the found k-mer in\nthe query and those in the unitig. For each matching character\nfound, local offset and k-mer end, i, are augmented by 1. This\nwalk ends when one of the following three conditions is met:\n1. the end of the query is reached; 2. the end of the unitig\nis reached; 3. a mismatch is found. In the latter two cases,\nprev is set back to −∞ and the search continues as described\nabove from the last k-mer not yet localized.\nIt is worth mentioning that in many applications, a k-mer\nis considered equal to its reverse complement. Thus if a k-\nmer Y is not found in the SBWT of the DSPSS it could still\nbe possible to find its reverse complement. We can take care\nof this by either doubling the query time by querying the k-\nmer in both orientations, or by up to doubling the index size\nby indexing the reverse complements of the unitigs as well.\nWe can store a permutation mapping the reverse unitigs to\ntheir forward counterparts, together with a bit vector marking\nthe starting k-mers of the reversed unitigs to easily check the\nunitig orientation. This allows us to report the locations of the\nfound k-mers in the forward unitigs.\nV. U NITIG ORIENTATION\nDue to the double strandedness of DNA, the existence of\na k-mer in a genome implies the existence of its reverse\ncomplement on the opposite strand. A common space-saving\ntrick is to index only one orientation of each k-mer, and con-\nsider the other orientation automatically present. For example,\nunitig extraction tools like BCALM2 [10], C UTTLEFISH 2\n[23] and GGCAT [11] by default produce only canonical\nunitigs, which are defined to be those unitigs which are\nlexicographically smaller or equal to their reverse complement.\nIn the best case, this reduces the number of k-mers to be\nindexed by a factor of 2. However, when using the SBWT to\nindex the k-mers in the unitigs, there is another complication\narising from the need to add incoming dummy node paths to\nall k-mers that do not have a predecessor in the de Bruijn\ngraph. If the unitigs are oriented in a canonical way as\ndescribed above, then due to chance we can expect adjacent\nunitigs will roughly half the time get oriented in opposite\ndirections, and so the number of unitigs (and therefore k-mers)\nwithout a predecessor may be unnecessarily large compared\nto a graph where neighbors are oriented in the same direction.\nHowever, since the de Bruijn graph can have cycles, it may\nbe unavoidable that at least some neighbors end up facing in\nopposite directions, and an interesting combinatorial optimiza-\ntion problem to orient the unitigs emerges.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n8\nFig. 3. A schematic representation of the proposed data structure for localization queries on the running example. Given the k-mer index\nwe compute for each position of the query, from left to right, the shortest unique suffixes array SFSk,1, which gives us the shortest unique\nfinimizers SUF, and the k-bounded matching statistics MSk. To the left we illustrate the ExtendRight operation, with extending character G\nand (magenta) initial interval [3, 5]AA, which leads to the (teal) interval [11, 11]AAG, further extended to [17, 17]AAGT(violet). We can perform\na single ContractLeft operation without modifying this interval. A further left contraction would lead to the (black) interval [17, 18]GT. The\nflow of the algorithm is represented with arrows. From MSk we can identify the rightmost branch for each k-mer in the query, if it exists.\nThe k-mer (orange, i = 4) AAGT is found in a non-branching path in the DBG, thus we have to look up the finimizer dictionary. The\nfinimizer is AGG and its colexicographic interval in the SBWT is [11, 11]. The global offset of the finimizer is obtained with a rank query\non F , G[F.rank(11)] = 5. We can get the unitig identifier with a predecessor query on E (unitig index), E.pred(5) = 1. The starting point\nof the k-mer is 3 = 5 − 0 + (4 − 3) − k + 1 Finally AAGT is localized in U ′ with result (1, 3). The localization query works differently\nfor GT AA(blue, i = 6) as we can identify a branch, at index 6, thus we know that the k-mer appears at the beginning of a unitig, which\nis identified with a rank query on U start, U start.rank(5) = 1. The starting point of the k-mer is 1 = 6 − 6 + 1. This last k-mer of the\nquery is then localized in U ′ with output (1, 1).\nIn this section, we analyze the space overhead caused by\nthe dummy nodes, and describe a heuristic to orient the unitigs\nto minimize the overhead. A similar heuristic is implemented\nin Metagraph [22] to reduce the number of dummy nodes in\nthe BOSS data structure. The optimal solution to the problem\nis left as an open problem.\nA. Dummy node overhead\nLet R be a set of k-mers and R′ ⊆ R be the subset of\nk-mers such that X ∈ R′ iff ̸ ∃Y ∈ R such that X[1..k −1] =\nY [2..k]. Let D = S\nx∈R′ {x[1..i] | i = 1 , . . . , k − 1} be the\nset of distinct proper prefixes of R′. Each element in D adds\na node to the SBWT graph (Figure 1, right). The exact space\ncost of this depends on the SBWT encoding used [3], [4]. For\nsimplicity here we assume each node has a unit cost, and edges\nare free. Ignoring lower-order terms in the space complexity of\nrank queries, this model is accurate for the bit matrix encoding\nof the SBWT used for experiments in this paper.\nThe number of dummy nodes introduced is equal to the\nnumber of nodes in the trie of the elements in D. Let x1, . . . xn\nbe the list of elements of R′ in lexicographic (note: not\ncolexicographic) order. Under the cost model defined above,\nthe total cost is |x1| − 1 + Pn\ni=2 (k − 1 − lcp(xi, xi−1))) =\nn(k − 1) − Pn\ni=2 lcp(xi, xi−1), where lcp(xi, xi−1) is the\nlength of the longest common prefix of k-mers xi and xi−1.\nB. Heuristic\nAs an exact solution seems difficult, we simplify the prob-\nlem by just trying to minimize n = |R′|. That is, minimizing\nthe number of unitigs without a predecessor. Our heuristic\nis best described as a series of breadth-first searches on the\nbidirected overlap graph of the unitigs. In a bidirected graph,\neach endpoint of an edge has an orientation, denoted with ⊕ or\n⊖. That is, edges are defined as ordered 4-tuples (u, v, s u, sv),\nwhere u and v are nodes and su and sv are the orientations\nof the endpoints. We use these kinds of edges to model all\nthe possible connections between unitigs and their reverse\ncomplements. Given a unitig v, we denote with v⊖ the reverse\ncomplement of v, and with v⊕ the unitig v itself, and we define\nthat there is an edge (u, v, s u, sv) iff the last k − 1 characters\nof usu are equal to the first k − 1 characters of vsv. Note that\nwith this definition, the existence of edge (u, v, s u, sv) implies\nthe existence of the symmetric edge (v, u, ¬sv, ¬su), where ¬\nflips the sign of an edge. We call a node u terminal iff for\nevery edge (u, v, su, sv), the sign su is the same.\nA breadth-first search maintains a queue of pairs (u, su),\nwhere u is a unitig and su is the orientation of the unitig. The\nqueue is initialized by pushing an arbitrary terminal node and\nassigning a forward orientation to it. Then, at each iteration,\nthe queue is popped, and if the node of the obtained pair\n(u, su) has not been visited yet with either orientation, we\nassign orientation su to u, and push all its neighbors to the\nqueue in orientations in which they overlap with (u, su). If the\nqueue becomes empty, but not all terminal nodes have been\nvisited so far, we push an arbitary unvisited terminal node with\na forward orientation and continue. If all terminal nodes have\nbeen visited, we push an arbitrary node. The algorithm con-\ncludes when all nodes have been visited. An implementation\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n9\nis available at https://github.com/jnalanko/unitig_flipper.\nVI. E XPERIMENTS\nExperimental Machine. Our experimental machine had four\nIntel Xeon E7-4830 v3 CPUs, each running at 2.10 GHz with\n12 cores (48 cores total); 30 MiB of L3 cache, and 1.5 TiB of\nRAM and a 12 TiB serial ATA hard disk. The OS was Linux\n(Ubuntu 18.04.5 LTS) kernel version 5.4.0-58-generic, with\ncompiler g++ version 10.3.0 (flags -O3 -march=native\nand -DNDEBUG). Runtimes were measured by instrumenting\ncode with calls to the high-resolution clock (std::chrono).\nDatasets. We experimented on three different datasets that\nrepresent typical targets for k-mer indexing (see also Table I).\n1) A pangenome of 3682 E. coli genomes 5.\n2) A set of 17,336,887 Illumina HiSeq 2500 reads of\nlength 502 sampled from the human gut (SRA identifier\nERR5035349) in a study on bile acid malabsorption\n[21].\n3) A set of 501,216 Nanopore sequencing reads of E. coli\ngenomes (SRA identifier SRR25689478) in a study on\nassorted bacteria with public health implications 6.\nA. Finimizer Length and Density\nFigure 4 shows the distribution of lengths of shortest unique\nfinimizers in our three datasets with k = 31. The distributions\nhave a peak around length 14 with a long tail to the right. The\npeak shifts to the right as the number of k-mers in the dataset\ngrows. Over 5% of the finimizers in the Nanopore dataset had\nthe longest possible length 31. We think this is due to the\nhigh error rate of Nanopore reads: If a true k-mer is read by\nthe sequencer at least three times, once without errors, once\nwith an error in the first position and once with an error just\nbefore the first position, then neither of its two (k − 1)-mers\nis DBG-unique, and the only possible finimizer is the k-mer\nitself.\nFigure 5 shows various statistics on shortest t-bounded\nfinimizers as a function of t in the 3682 E. coli genomes\ndataset. We find that while the number of distinct finimizers\ndecreases with increasing values of t, the total number of DBG\noccurrences grows. This means that the choice t = 1 provides\nboth a smaller index size (since we store a single integer\nper occurrence) and faster queries as in the worst case the\nquery needs to check all occurrences of a finimizer. However,\nthe average finimizer length decreases when t grows, which\nmeans that indexing strategies whose space depends on the\nlengths of finimizers (such as those based on the Aho-Corasick\nautomaton) might favor larger choices of t.\nA fixed-length minimizer scheme that orders m-mers ran-\ndomly has a probability of 2/(w + 1) of selecting a position\nas the ending position of a minimizer, where w = k −m+1 is\nthe number of m-mers inside a k-mer [37]. Experimentally, we\nobserved that on random data, this formula also approximately\nholds for finimizers so that w is equal to k − ¯m + 1, where ¯m\nis the average length of chosen finimizers on the dataset: On a\n5The collection is available at zenodo.org/record/6577997.\n6See https://www.ncbi.nlm.nih.gov/bioproject/812595.\n0 5 10 15 20 25 30\nFinimizer length\n0.00\n0.05\n0.10\n0.15\n0.20\n0.25\n0.30\n0.35\n0.40\nFrequency\nFinimizer length distribution\nMetagenome short reads\nE. coli Nanopore reads\n3682 E. coli genomes\nFig. 4. Length distribution of the distinct shortest unique finimizers on three\ndifferent datasets with k = 31.\nrandom DNA string consisting of 106 characters, with k = 63,\nthe average finimizer length was 9.50, and the shortest unique\nfinimizer scheme picked 3.76% positions as ending positions\nof finimizers, whereas the formula predicts a density of 3.60%.\nHowever, on non-random data, the gap becomes wider. On the\nunitigs of a single E. coli genome (average unitig length 6048)\nfrom our 3682 E. coli genomes dataset, 4.3% of positions were\nselected as ending positions of finimizers, but with average\nlength ¯m = 10.50, the formula predicts a density of only\n3.67%.\nB. Index Size and Query Performance\nWe have developed a prototype k-mer localization index\nbased on Shortest Unique Finimizers, called F INITO . We now\nreport on performance benchmarks of this prototype on three\nlarge datasets. We use SSH ASH [33] v3.0.0 (commit 220b2c0)\nas our performance baseline because it was recently shown to\nclearly outperform other fast and popular k-mer localization\ntools (e.g., BL IGHT [28])—it is the fastest k-mer localization\nmethod we know 7.\nConstruction. The input was preprocessed into canonical\nunitigs before index construction, using GGCAT [11]. The\ntime for this preprocessing in not included in the construction\ntimes. We then employed a heuristic method to reduce dummy\nnodes in the SBWT on all datasets. The method 8 is described\nin the supplementary materials, and may be of independent\ninterest.\nFrom here, for k = 31 and m = 20 SSH ASH ,with\notherwise default parameters, took 11, 3, and 2 minutes to\nconstruct, respectively, for the Metagenome, Nanopore, and\nE.Coli datasets, with peak memory of 11GB, 4GB, and 2GB.\nFor the same k, construction of the double-stranded F INITO\nindex was longer, resp.: 299, 58, and 13 minutes (peak\nmemory, resp.: 140GB, 42GB, and 12GB), more than half\n7As described in [33], SSH ASH only claims to solve the problem of k-mer\nlookup, but the implementation requires only trival modifications to report\nk-mer localization.\n8See https://github.com/jnalanko/unitig_flipper\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n10\n0 2 4 6 8 10\nt\n0\n1\n2\n3# distinct finimizers\n×107 # Distinct finimizers\n0 2 4 6 8 10\nt\n0\n10\n20\n30Average length\nAverage finimizer length\n0 2 4 6 8 10\nt\n0.0\n0.1\n0.2\n0.3\n0.4\n0.5\nSampling density\nDBG sampling density\n0 2 4 6 8 10\nt\n0\n1\n2\n3\n4\n5Average frequency\nAverage finimizer frequency\nk = 21\nk = 31\nk = 63\nk = 127\nFig. 5. Statistics on shortest t-bounded finimizers on the 3682 E. coli genomes dataset. The DBG sampling density at bottom left is defined as the number\nof DBG occurrences of all distinct finimizers, divided by the number of distinct nodes in the DBG. The frequency plotted at the bottom right is the average\nnumber of DBG occurrences of finimizers.\nTABLE I\nStatistics on the raw genomic datasets.k-mer counts considers canonicalk-mers. We derived unitigs for each dataset and built a SBWT. $\noriginal is the number of dummy nodes in the original unitig ordering, and $ flipped is the number after applying the unitig orientation\nheuristic described in Section V.\nSequences Total length Unique 31-mers Unitigs Unitig length $ original $ flipped\nE. coli 745,409 18,957,578,183 170,648,610 8,939,732 438,840,570 52,528,157 4,838,079\nMetagenome 17,336,887 8,703,117,274 2,761,523,935 76,958,596 4,182,332,581 517,193,809 191,843,337\nNanopore 501,216 2,888,677,990 483,259,051 29,581,392 1,370,700,811 119,248,987 10,994,322\nof which was spent on SBWT construction9. We remark that\nalthough F INITO construction time is higher than SSH ASH, it\nstill acceptably fast in real terms. Moreover, the construction of\nFINITO currently runs on a single thread, whereas SSH ASH\nuses 8. The bottleneck in our construction process—SBWT\ncomputation—could be significantly sped up with paralleliza-\ntion, a task we leave for future work.\nQueries. We generated two types of query sets: random\nsequences of length 200 (negative queries ), the k-mers of\nwhich are highly unlikely to have any matches in the datasets;\nand positive queries were generated by randomly sampling\n200-mers from each dataset. We measured the time F INITO\nand SSH ASH took to look up every k-mer in 500,000 positive\nand negative length-200 sequences, for various values ofk (see\nFigure 6). To match Sshash queries, we consider ak-mer found\nif it is found in either forward or reverse complement orienta-\ntion. As mentioned before, we achieved this either by using an\nindex that contains both the forward and reverse complement\nstrains of all unitigs, and storing the permutation mapping\nreverse and forward counterparts to report the localizations\n9Construction time for the single-stranded index were faster, resp. 131, 26,\n8 minutes (using resp. 4.6, 2, 2.5GB of RAM), with a corresponding slow\ndown in query time.\non the original DNA strand, or by using an index containing\nonly one direction, but running all queries in both directions.\nThe former approach being roughly twice as fast but taking\ntwice the space as the latter. The index size of F INITO was\ntypically larger than Sshash, but the query times of both\ntools were similar, within a small constant factor from each\nother, depending on the configuration and the dataset, with\nnegative queries favoring F INITO . Alanko et al. [4] provide\nseveral data structures to support ExtendRight queries on the\nSBWT. We use their so-called Plain Matrix representation.\nOther data structures could reduce the size of our index at the\nexpense of an increased query time, at least with current SBWT\nincarnations. In this respect, it is worth noting that SBWT data\nstructures are still very new, and that further engineering may\nyet lead to significant performance gains.\nVII. C ONCLUSIONS AND FUTURE WORK\nIn the 20 years since minimizers emerged, researchers have\ninvented various heuristics to avoid very frequent minimizers.\nIn this paper we have demonstrated an alternative solution\nto the problem that works by letting the minimizer length\nvary. This way, when indexing DSPSS representations, we\ncan get a guarantee for the maximum frequency and thus\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n11\n0 20 40 60\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\nTime (µs/kmer)\n12\n16 20\n12\n1620\nE. coli\n0 20 40 60\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\n12\n16\n20\n12\n16\n20\nNanopore\n0 20 40 60\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\n12\n16\n12\n16\n20\nMetagenome\nSshash\nFinito\nFinito+rev\nk=21\nk=31\n0 20 40 60\nIndex size (bits/kmer)\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\nTime (µs/kmer)\n12\n16\n20\n1620\n0 20 40 60\nIndex size (bits/kmer)\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\n16\n20\n1620\n0 20 40 60\nIndex size (bits/kmer)\n0.0\n0.5\n1.0\n1.5\n2.0\n2.5\n3.0\n3.5\n4.0\n16\n12\n16\n20\nFig.\n6. Query performance of SSH ASH and F INITO . Top: Time and index size for positive queries for varying k and m on each dataset.\nBottom: Time and index size for negative queries.SShash results for Metagenome with k = 21, m = 20 are not plotted because construction\nfailed due to collisions of 64-bit hash values.\nmaximum localization query time. A distinct advantage the\nShortest Unique Finimizers, we have introduced, have over\nminimizers is that they are parameter free.\nWhile our current prototype is already competitive with\nstate-of-the-art minimizer-based localization methods, we be-\nlieve it can be optimized further, e.g., by adding word-\nlevel parallelism to our ContractLeft implementation, and by\nreducing the time and space for ExtendRight by engineering\nthe structures from [4] to use faster recent predecessor [13]\nand rank [3], [8] structures.\nThe SBWT offers a convenient way to index finimizers for\nlookup, however more direct methods based on, for example,\nmultiple pattern matching data structures (e.g., [39]), may\nlead to smaller and faster indexes. The principle challenge\n(compared to minimizers) is to deal with variable-length keys.\nIn this respect, exploring trie-based solutions may be fruitful.\nFinally, interesting applications for finimizers beyond k-\nmer localization include: skipping heuristics in pseudoalign-\nment [6]; approximate k-mer matching; and replacing hashing\nwith finimizers in spectrum-preserving tilings [16].\nA C OMPETING INTERESTS\nNo competing interest is declared.\nB A UTHOR CONTRIBUTIONS STATEMENT\nAll authors contributed equally.\nC A CKNOWLEDGMENTS\nThis work was supported in part by the Academy of Finland\nvia grants 339070 and 351150.\nREFERENCES\n[1] P. Abedin, M. O. Külekci, and S. V . Thankachan. A survey on shortest\nunique substring queries. Algorithms, 13(9):224, 2020.\n[2] J. N. Alanko, E. Biagi, and S. J. Puglisi. Longest common prefix arrays\nfor succinct k-spectra. In Proc. SPIRE, LNCS 14240, pages 1–13.\nSpringer, 2023.\n[3] J. N. Alanko, E. Biagi, S. J. Puglisi, and J. Vuohtoniemi. Subset wavelet\ntrees. In Proc. SEA, LIPIcs 265, pages 4:1–4:14. Schloss Dagstuhl, 2023.\n[4] J. N. Alanko, S. J. Puglisi, and J. Vuohtoniemi. Small searchable\nk-spectra via subset rank queries on the spectral Burrows-Wheeler\ntransform. In Proc. ACDA, pages 225–236. SIAM, 2023.\n[5] D. Belazzougui and F. Cunial. Indexed matching statistics and shortest\nunique substrings. In Proc. SPIRE 2014, pages 179–190. Springer, 2014.\n[6] N. L. Bray, H. Pimentel, P. Melsted, and L. Pachter. Near-optimal\nprobabilistic RNA-seq quantification. Nature Biotechnology, 34(5):525–\n527, 2016.\n[7] R. Cánovas and G. Navarro. Practical compressed suffix trees. In\nP. Festa, editor, Proc. 9th International Symposium Experimental Al-\ngorithms (SEA), volume 6049 of Lecture Notes in Computer Science,\npages 94–105. Springer, 2010.\n[8] M. Ceregini, F. Kurpicz, and R. Venturini. Faster wavelet trees with\nquad vectors. CoRR, abs/2302.09239, 2023.\n[9] R. Chikhi, A. Limasset, S. Jackman, J. T. Simpson, and P. Medvedev.\nOn the representation of de Bruijn graphs. In Proc. RECOMB, LNCS\n8394, pages 35–55. Springer, 2014.\n[10] R. Chikhi, A. Limasset, and P. Medvedev. Compacting de Bruijn graphs\nfrom sequencing data quickly and in low memory. Bioinformatics,\n32(12):i201–i208, 06 2016.\n[11] A. Cracco and A. Tomescu. Extremely fast construction and querying\nof compacted and colored de Bruijn graphs with GGCAT. Genome res.,\n05 2023.\n[12] S. Deorowicz, M. Kokot, S. Grabowski, and A. Debudaj-Grabysz. KMC\n2: fast and resource-frugal k-mer counting. Bioinformatics, 31(10):1569–\n1576, 2015.\n[13] D. Díaz-Domínguez, S. Dönges, S. J. Puglisi, and L. Salmela. Simple\nruns-bounded FM-index designs are fast. In Proc. SEA, LIPIcs 265,\npages 7:1–7:16. Schloss Dagstuhl, 2023.\n[14] B. Ekim, B. Berger, and R. Chikhi. Minimizer-space de bruijn graphs:\nWhole-genome assembly of long reads in minutes on a personal com-\nputer. Cell Systems, 12(10):958–968.e6, 2021.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n12\n[15] M. Erbert, S. Rechner, and M. Müller-Hannemann. Gerbil: a fast\nand memory-efficient k-mer counter with GPU-support. Algorithms for\nMolecular Biology, 12(9), 2017.\n[16] J. Fan, J. Khan, G. E. Pibiri, and R. Patro. Spectrum preserving tilings\nenable sparse and modular reference indexing. In Proc. RECOMB,\nLNCS 13976, pages 21–40. Springer, 2023.\n[17] U. Ferraro-Petrillo, M. Sorella, G. Cattaneo, R. Giancarlo, and S. E.\nRombo. Analyzing big datasets of genomic sequences: fast and scalable\ncollection of k-mer statistics. BMC Bioinformatics, 20(Suppl 4):138,\n2019.\n[18] J. Fischer. Combined data structure for previous-and next-smaller-values.\nTheoretical Computer Science, 412(22):2451–2456, 2011.\n[19] G. Holley and P. Melsted. Bifrost: highly parallel construction and\nindexing of colored and compacted de Bruijn graphs. Genome biology,\n21(1):1–20, 2020.\n[20] C. Jain et al. Weighted minimizer sampling improves long read mapping.\nBioinf., 36(Suppl 1):i111–i118, 07 2020.\n[21] I. B. Jeffery et al. Differences in fecal microbiomes and metabolomes\nof people with vs without irritable bowel syndrome and bile acid\nmalabsorption. Gastroenterology, 158(4):1016–1028, 2020.\n[22] M. Karasikov, H. Mustafa, D. Danciu, M. Zimmermann, C. Barber,\nG. Rätsch, and A. Kahles. Metagraph: Indexing and analysing nucleotide\narchives at petabase-scale. BioRxiv, 2020.\n[23] J. Khan, M. Kokot, S. Deorowicz, and R. Patro. Scalable, ultra-fast, and\nlow-memory construction of compacted de Bruijn graphs with Cuttlefish\n2. Genome Biology, 23, 09 2022.\n[24] H. Li. Aligning sequence reads, clone sequences and assembly contigs\nwith bwa-mem. arXiv preprint arXiv:1303.3997, 2013.\n[25] H. Li. Minimap2: pairwise alignment for nucleotide sequences. Bioin-\nformatics, 34(18):3094–3100, 05 2018.\n[26] C. Marchet, C. Boucher, S. J. Puglisi, P. Medvedev, M. Salson, and\nR. Chikhi. Data structures based on k-mers for querying large collections\nof sequencing data sets. Genome Res., 31(1):1–12, 2021.\n[27] C. Marchet, Z. Iqbal, D. Gautheret, M. Salson, and R. Chikhi. Reindeer:\nefficient indexing of k-mer presence and abundance in sequencing\ndatasets. Bioinf., 36:i177–i185, 07 2020.\n[28] C. Marchet, M. Kerbiriou, and A. Limasset. BLight: efficient exact\nassociative structure for k-mers. Bioinformatics, 37(18):2858–2865,\n2021.\n[29] G. Navarro. Compact Data Structures – A practical approach. Cam-\nbridge University Press, 2016.\n[30] J. Nyström-Persson, G. Keeble-Gagnère, and N. Zawad. Compact and\nevenly distributed k-mer binning for genomic sequences. Bioinformatics,\n37(17):2563–2569, 2021.\n[31] E. Ohlebusch, S. Gog, and A. Kügel. Computing matching statistics\nand maximal exact matches on compressed full-text indexes. In Proc.\nSPIRE 2010, pages 347–358. Springer, 2010.\n[32] Y . Orenstein et al. Designing small universal k-mer hitting sets for\nimproved analysis of high-throughput sequencing. PLoS Computational\nBiology, 13(10):e1005777, 2017.\n[33] G. E. Pibiri. Sparse and skew hashing of K-mers. Bioinformatics,\n38(Suppl 1):i185–i194, 06 2022.\n[34] A. Rahman and P. Medevedev. Representation of k-mer sets using\nspectrum-preserving string sets. J. Computational Biology, 28(4):381–\n394, 2021.\n[35] M. Roberts et al. A preprocessor for shotgun assembly of large genomes.\nJ. Computational Biology, 11(4):734–752, 2004.\n[36] M. Roberts, W. Hayes, B. R. Hunt, S. M. Mount, and J. A. Yorke.\nReducing storage requirements for biological sequence comparison.\nBioinformatics, 20(18):3363–3369, 07 2004.\n[37] S. Schleimer, D. S. Wilkerson, and A. Aiken. Winnowing: local\nalgorithms for document fingerprinting. In Proceedings of the 2003\nACM SIGMOD international conference on Management of data, pages\n76–85, 2003.\n[38] S. Schmidt, S. Khan, J. N. Alanko, G. E. Pibiri, and A. I. Tomescu.\nMatchtigs: minimum plain text representation of k-mer sets. Genome\nBiology, 24(1):1–32, 2023.\n[39] S. Wu and U. Manber. A fast algorithm for multi-pattern searching.\nTechnical Report TR94-17, University of Arizona. Department of Com-\nputer Science Tucson, AZ, 1994.\n[40] H. Zheng, G. Marçais, and C. Kingsford. Creating and using minimizer\nsketches in computational genomics. J. Computational Biology, 30(0):1–\n26, 2023.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\nSupplementary material for Finimizers: Variable-length\nbounded-frequency minimizers for k-mer sets\nJarno N. Alanko, Elena Biagi, Simon J. Puglisi\n1 k-bounded Matching Statistics\nAlgorithm 1 k-bounded matching statistics.\nInput: A query string Q and a data structure supporting ExtendRight and ContractLeft on the padded k-mer spectrum.\nOutput: MSk[1..|Q|]\n1: [ℓ, r] ← [1, n] ▷ The current colexicographic interval\n2: d ← 0 ▷ The length of the current match, up to k\n3: MSk[1..|Q|] ← Array of length |Q|.\n4: for i = 1..|Q| do\n5: while d > 0 and ExtendRight([ℓ, r], Q[i]) = ∅ do\n6: [ℓ, r] ← ContractLeft([ℓ, r])\n7: d ← d − 1\n8: if ExtendRight([ℓ, r], Q[i]) ̸= ∅ then\n9: [ℓ, r] ← ExtendRight([ℓ, r], Q[i])\n10: d ← min(k, d + 1)\n11: MSk[i] ← (d, r − ℓ + 1, ℓ) ▷ Len, Freq, Colex start\nDefinition 1. (k-bounded matching statistics) The k-bounded matching statistics of a query sequence Q[1..m], versus\na padded k-spectrum R+, is the array MS k[1..m] such that MS k[i] = (di, pi, ℓi), where di is the largest non-negative\ninteger such that Q(i − di..i] is a suffix of at least one k-mer in R+, pi is the DBG frequency of Q(i − di..i] and ℓi\nis the start of the colexicographic interval of Q(i − di..i].\nTheorem 1. Given constant-time ContractLeft and ExtendRight queries on the padded k-spectrum R+, Algorithm 1\ncomputes M Sk for query Q in time O(|Q|).\nProof. The algorithm maintains the invariant that at the end of iteration i, the interval [ℓ, r] is the interval of the\nlongest suffix of Q[1..i] that occurs as suffix of at least one k-mer. If [ ℓ, r] is the longest match at the end of iteration\ni − 1, then at iteration i after lines 5 – 7, we have the interval of the longest match ending at i − 1 that can be\nextended with Q[i], or the interval of the empty string if Q[i] does not exist in any k-mer. If the longest match\nending at i is nonempty, it is a right extension of such a string, so the invariant is maintained. The algorithm works\nin O(|Q|) time because d is always non-negative, and the while-loop can run only d times because d is incremented\nat most |Q| times, and every run of the while-loop decrements d.\n1\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n2 Shortest Frequency-bounded Suffixes\nAlgorithm 2 Shortest frequency-bounded suffix statistics.\nInput: A query string Q, frequency upper bound t, and a data structure supporting ExtendRight and ContractLeft on the\npadded k-mer spectrum.\nOutput: SFSk,t[1..|Q|]\nAssumes all characters of the query Q exist in the index at least once.\n1: [ℓ, r] ← [1, n] ▷ The current colexicographic interval\n2: d ← 0 ▷ The length of the current match, up to k\n3: SFSk,t[1..|Q|] ← Array of length |Q| initialized to values ⊥\n4: for i = 1..|Q| do\n5: while ExtendRight([ℓ, r], Q[i]) = ∅ do\n6: [ℓ, r] ← ContractLeft([ℓ, r])\n7: d ← d − 1\n8: [ℓ, r] ← ExtendRight([ℓ, r], Q[i])\n9: d ← min(k, d + 1)\n10: while |ContractLeft([ℓ, r], 1)| ≤ t do\n11: [ℓ, r] ← ContractLeft([ℓ, r])\n12: d ← d − 1\n13: if r − ℓ + 1 ≤ t then\n14: SFSk,t[i] ← (d, r − ℓ + 1, ℓ)\nDefinition 2. (Shortest frequency-bounded suffixes) The shortest suffix statistics of a query sequenceQ[1..m], versus\na padded k-spectrum R+, is the array SF Sk,t[1..m] such that SF Sk,t[i] = (di, pi, ℓi), where di is the smallest positive\ninteger such that Q(i − di..i] is a suffix of pi (pi ≤ t) k-mers in R+, and ℓi is the start of the colexicographic interval\nof Q(i − di..i]. If no such di exists, SF Sk,t[i] = ⊥.\nTheorem 2. Given constant-time ContractLeft and ExtendRight queries on the padded k-spectrum R+, Algorithm 2\ncomputes SF Sk for query Q in time O(|Q|).\nProof. The algorithm maintains the invariant that at the end of iteration i, the interval [ℓ, r] is the interval of the\nshortest suffix of Q[1..i] that occurs as suffix of at most t k -mers, if exists, or otherwise [ℓ, r ] is the longest suffix of\nQ[1..i] that occurs in the k-mers. To see this, we consider two cases.\n• (i) If [ ℓ, r] is the interval of the longest match at the end of iteration i − 1, then in iteration i, the lines 5 –\n9 update it to the longest match up to position i. If the interval becomes of size ≤ t after lines 5 – 9, then\nlines 10 – 12 contract it to the shortest match with frequency at most t and the invariant continues to hold.\nOtherwise if the interval size did not shrink to at most t after lines 5 – 9, then lines 10 – 12 are not run and\nthe interval remains the longest match. The invariant continues to hold.\n• (ii) If [ ℓ, r] is the shortest match with frequency at most t at the end of iteration i − 1, then\n– (iia) if the loop on lines 5 – 7 is run at least once, it becomes the longest match after the right extension\nat line 8 (a left contraction was needed for a right extension, so the positions that were contracted can\nnot be a part of the longest match). If the loop on lines 10 – 12 is run at least once, it is contracted\nto the shortest match of frequency at most t, otherwise the interval is not changed and it remains the\nlongest match, maintaining the invariant.\n– (iib) if the loop on lines 5 – 7 is not run at least once, then the interval remains unchanged up to line 8,\nand the right extension on line 8 can only make the frequency smaller, so the frequency remains at most\nt. Lines 10 – 12 contract it to the shortest match with frequency at most t, maintaining the invariant.\nIf the assignment on line 14 happens, it is correct since then [ l, r] is of size at most t and by lines 10 – 12 the length\nd is minimal. If the assignment does not happen, then SF Sk[i] is correctly left as ⊥ since now the interval [ ℓ, r]\nrepresents the longest available match ending at i, and thus a match with frequency t or less can not exist, as it\nwould have to be longer than the longest match.\nThe algorithm works in O(|Q|) time because the while-loop conditions are false at most |Q| times each, and\nevery time they evaluate to true, the left contraction inside can be charged to a distinct right extension on line 8,\nwhich is executed |Q| times.\n2\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint \n\n3 Localization Proof of Correctness\nSuppose we have query a Q such that k-mer α = Q(i − k..i] is known to exist in the DBG according to the\nmatching statistics information. We now prove that our algorithm localizes α to the correct unitig.\n• Case 1: Localization via the branch dictionary . Suppose there exists a substring β = Q[p..q] with\ni − k + 1 ≤ p ≤ q ≤ i such that the DBG frequency of β is 1, and β is the suffix of a k-mer that is at the\nstart of a unitig X. Further, suppose that q is the largest possible (up to i) with this property. Then, our\nalgorithm localizes α to X. Let δ = Q[q + 1..i]. First, we claim that βδ is a substring of X.\n– Proof. For a contradiction, suppose this is not the case. Then, since β is unique in the DBG and α is\nknown to exist in the DBG, it must be the case that X ends before the end of δ, that is, |X| < k + |δ|,\nwhich means that the path labeled with βδ contains an edge (u, v) such that u has outdegree larger than\n1. Let βδ[1..r] be the path label up to and including v. This path label is still unique in DBG, and it is\na suffix of the first k-mer of its unitig. But then q was not the largest possible since position q + r also\nfulfills the conditions imposed on q, contradiction.\nLet α = γβδ . We claim that γβ is a substring of X.\n– Proof. Since α exists in the DBG and β is unique in the DBG, the occurrence of β in the graph has\nan incoming path with label γ. Since β = X[k − |β| + 1..k], we have that k − |β| − |γ| + 1 ≥ 1 and so\nX[k − |β| − |γ | + 1..k] = γβ.\nSince β is DBG-unique, and both γβ and βδ are in X, we conclude that α = γβδ must also be in X.\n• Case 2: Localization via the finimizer dictionary. Suppose there does not exist β in the case above.\nLet ρ be the shortest unique finimizer of α, so α = τ ρτ′. Then our algorithm localizes α to the unitig which\ncontains an occurrence of ρ the furthest from the start of the unitig. First, we claim that the path of α is\nnon-branching DBG.\n– Proof. Let v1, . . . , vk+1 be the path labeled with α. We consider incoming branches and outgoing\nbranches as separate cases:\n∗ Suppose there exists j > 1 such that the indegree of vj is at least 2. Then, ρ is the suffix of two\ndistinct k-mers, so it is not unique in the DBG, a contradiction.\n∗ Suppose there exists j < k + 1 such that the outdegree of vj is at least 2. Then, vj must be after\nthe ending point of ρ on the path, since otherwise ρ is not DBG-unique. But then the label of the\nsubpath from the start of ρ to the node vj+1 satisfies the conditions for β, which is a contradiction\nsince such β was assumed to not exist.\nThe global offset stored for ρ is the one furthest from the unitig start. We claim that the global offset of ρ\nthat we have stored is in the unitig that contains α:\n– Proof. The string ρτ ′ is fully contained in the same unitig as ρ, because otherwise the path of ρτ ′ would\nbe branching, but that cannot be that case since α = τ ρτ′ in non-branching. Since we are at the unitig\nX that contains ρ the furthest to the right as possible, the whole string α = τ ρτ′ must be found fully in\nX.\n3\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted February 21, 2024. ; https://doi.org/10.1101/2024.02.19.580943doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}