Martin C Frith

Laurent Noé

Gregory Kucherov

Minimally-overlapping words for sequence similarity search

Analysis of genetic sequences is usually based on finding similar parts of sequences, e.g. DNA reads and/or genomes. For big data, this is typically done via "seeds": simple similarities (e.g. exact matches) that can be found quickly. For huge data, sparse seeding is useful, where we only consider seeds at a subset of positions in a sequence.

Here we study a simple sparse-seeding method: using seeds at positions of certain "words" (e.g. ac, at, gc, or gt). Sensitivity is maximized by using words with minimal overlaps. That is because, in a random sequence, minimally-overlapping words are anti-clumped. We provide evidence that this is often superior to acclaimed "minimizer" sparseseeding methods. Our approach can be unified with design of inexact (spaced and subset) seeds, further boosting sensitivity. Thus, we present a promising approach to sequence similarity search, with open questions on how to optimize it.

Introduction 1.Seeds

Finding similar sequences, in large data, is typically done via "seeds": simple similarities that can be found quickly. The simplest type of seed is exact matches of a given length, e.g. 10 letters for DNA. The seed length affects the sensitivity and run time: shorter seeds are more sensitive, but find more hits that must then be checked. By lengthening the seeds, we can arbitrarily reduce the run time of the downstream steps, but not the time and memory usage for finding the seeds.

Sparse seeds

An alternative way to reduce time and/or memory use is sparse seeding. The simplest way is to only use seeds starting at every nth position in one of the two sequences being compared.

Note that, if we only use seeds at every nth position in both sequences, the sensitivity will be poor. E.g. if there are long similar segments, without insertions or deletions, starting at coordinate x in the first sequence and coordinate x + 1 in the second sequence, they are never hit if n > 1.

Sparse seeding reduces sensitivity, but we could then increase the sensitivity by shortening the seeds. This raises the prospect of reducing run time and/or memory use without loss of sensitivity.

Sparsity via words

An intriguing idea is to achieve sparsity by selecting seeds starting at positions of certain words. For example, if we only use seeds starting with a (Paul Horton, personal communication), we achieve 4-fold sparsity in both sequences without huge loss of sensitivity. We can imagine more complex variants, e.g. use seeds starting with any of these words: ac, at, gc, gt. Surprisingly, it makes a difference: it is better to use words that have minimal chance 1 of overlapping. That is because, in a random sequence, minimally-overlapping words occur with more uniform spacing, i.e. they are anti-clumped or under-dispersed; equivalently, their number of occurrences has lower variance.

Minimizers

A related idea is minimizers [SWA03, RHH + 04]. This method only uses seeds starting at positions that are "minimum" in any window of w consecutive positions (e.g. w = 7). Various definitions of "minimum" can be used; the simplest is: the sequence (i.e. suffix) starting at this position is alphabetically minimum. This is somewhat like using seeds starting with a. More complex orderings can be used, e.g. compare two suffixes using order c<a<t<g at oddnumbered bases and g<t<a<c at even-numbered bases, so that cgcg. . . is the minimum possible suffix [RHH + 04].

The resulting degree of sparsity is not obvious, and it depends on the ordering [MPB + 17]. Typically, a fraction 2/(w + 1) of positions is selected [START_REF] Schleimer | Winnowing: local algorithms for document fingerprinting[END_REF].

Another related idea is universal k-mer hitting sets [OPM + 17]. This means a set of length-k words, such that every possible length-L sequence contains at least one of the words. Recent studies have defined minimizer orderings based on universal kmer hitting sets, resulting in high sparsity for a given w [OPM + 17, MPB + 17, MDK18].

Minimizers have been described as "a central recent paradigm" [OPM + 17]: they have been widely used for sparse seeding (e.g. [Li18, JKD + 18]) and other applications (e.g. [LKH + 13, WS14, DKGDG15]).

Spaced and subset seeds

So far we have considered exact-match seeds, but inexact seeds are also used. One variant is spaced seeds, which allow mismatches at some fixed positions in the seed (e.g. positions 3 and 5 out of 9). Spaced seeds are often superior to exact seeds [START_REF] Ma | PatternHunter: faster and more sensitive homology search[END_REF], because their hits are less concentrated in overlapping clumps. Thus, spaced seeds have been designed by minimizing their "overlap complexity" [START_REF] Ilie | Multiple spaced seeds for homology search[END_REF], which is similar to minimizing the variance in number of hits [HLO + 16].

Subset seeds are a further generalization: they allow some mismatches (e.g. a↔g and c↔t) at fixed positions [START_REF] Noé | Improved hit criteria for DNA local alignment[END_REF]. This is useful for DNA, because a↔g and c↔t substitutions (termed "transitions") are often more frequent than the other types of substitution ("transversions"). Transition seeds have also been designed for use with every-nth sparsity [START_REF] Frith | Improved search heuristics find 20 000 new alignments between human and mouse genomes[END_REF].

Repeats

Natural DNA has many repeats, which are the main difficulty for similarity search. For example, a primate genome may have a million Alu elements, so naive comparison of such genomes yields an unmanageable 10 12 significant similarities. Our practical aim cannot be to find all significant similarities, but rather orthologs and/or strongest similarities. In any case, a seeding method must avoid getting too many repetitive seeds. One solution is to omit highfrequency seeds, another is to use variable-length seeds that are made longer until they are sufficiently rare [Csű04, KWS + 11].

Non-overlapping words

Since we are interested in minimally-overlapping words, let us consider non-overlapping words. A basic question is: what is the maximum possible number of non-overlapping words of some length k? That is, given an alphabet of size a (so there are a k possible words), what is the maximum possible number of words where no proper prefix of any word equals a proper suffix of any word? This seems hard to answer in general [START_REF] Simon | Non-overlapping codes[END_REF].

The following construction has been suggested for getting a large number of non-overlapping words [START_REF] Simon | Non-overlapping codes[END_REF]. Divide the alphabet into two subsets, e.g. {a} and {c, g, t}, and choose a prefix length j (0 < j < k). These words have no overlaps: words whose first j letters are from the 1st subset, whose (j +1)th and kth letters are from the 2nd subset, and whose letters between j + 1 and k have no run of ≥ j letters from the 1st subset.

Methods

The methods are available at: https://gitlab. com/mcfrith/noverlap

Mean and variance

Given a set of length-k words, let us consider their occurrence in a random i.i.d. length-s sequence. The expected number of total occurrences is:

E[X] = (s -k + 1)p , ( 1
)
where p is the total probability of any of the words occurring at a given position. The variance in occurrence number is:

Var[X] = E[X 2 ] -E[X] 2 . ( 2
)
To calculate this, let us define index variables I j as 1 if any of the words occurs at position j, else 0. So:

E[X 2 ] = E[(I 1 + I 2 + • • • + I s-k+1 ) 2 ] (3) = i,j E[I i I j ] . (4)
If we define l = |i -j|, then

E[I i I j ] =      p if l = 0 V,W B k-l V W P l V P k W if 0 < l < k p 2 if l ≥ k , (5)
where B m V W is defined to be 1 if the length-m suffix of word V equals the length-m prefix of word W , else 0. Also, P n W is the product of probabilities of the first n letters in word W . Thus, assuming that s ≥ 2k -2 (see the Supplement):

Var[X] = (s -k + 1)p + 2 V,W k-1 l=1 [(s -k + 1 -l)B k-l V W P l V P k W ] -[(2k -1)s -(3k -1)(k -1)]p 2 (6)
For circular sequences, the formulas are simpler (assuming s ≥ 2k -1):

E[X] = sp (7) Var[X] = s p -(2k -1)p 2 + 2 V,W k-1 l=1 (B k-l V W P l V P k W ) (8)
These formulas also apply to linear sequences when s k. With these formulas, the variance-to-mean ratio, also called index of dispersion, is independent of the sequence length. The formulas also simplify for linear sequences with s = 2k -1 (the smallest s where all kinds of pairwise overlap contribute):

E[X] = kp (9) Var[X] = kp -k 2 p 2 + 2 V,W k-1 l=1 [(k -l)B k-l V W P l V P k W ] (10)

Simulated sequences

To test homology detection, DNA sequences were simulated with the T92 model of evolution [START_REF] Tamura | Estimation of the number of nucleotide substitutions when there are strong transitiontransversion and G+C-content biases[END_REF]. This model has three input parameters: gc-content, transition/transversion rate ratio κ, and PAM substitution distance (Table 1).

For each test, 100 000 pairs of DNA sequences were simulated. The default parameters, unless specified otherwise, are: %g+c = 50, κ = 1 (unbiased), PAM = 20, sequence length = 100. A seeding method was deemed to find a pair of sequences if it found at least one match at identical coordinates of the pair.

To test specificity, two unrelated length-10 6 sequences were generated, and the number of seed pair matches counted. This is a proxy for the computational cost of checking all the seed hits.

Results

Non-overlapping DNA words

The maximum possible number of non-overlapping DNA words, for word length k = 2-6 (Table 2), was found by brute force clique search [KJ07]. For k < 6, Blackburn's construction [START_REF] Simon | Non-overlapping codes[END_REF] achieves this maximum. For k = 2, a maximum set is ry (r = a or g, y = c or t). In general, abb. . . (b = any base except a) is a good way to get non-overlapping words, and a nice generalization of Horton's idea. number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q 5 6 number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q PAM 50 q every-1st every-2nd every-4th every-8th every-16th

Every nth sparsity

We first tested every-nth sparsity (only using seeds starting at every nth position in one of the two sequences being compared), with exact-match seeds.

We defined "sensitivity" as % of sequence pairs with ≥ 1 seed match at homologous positions. As expected, if we increase sparsity without changing the seed length, both sensitivity and random hit count decrease (Figure 1). If we then shorten the seeds, the sensitivity and random hit count increase. The important result is that higher sparsity has lower sensitivity for a given random hit count. The exception is n = 2, which is no worse than n = 1, indeed giving us something for nothing: sparsity at no cost. A plausible explanation for why n = 1 is not better than n = 2 is that highly-overlapping seeds provide little independent information. This is also why spaced seeds are better than exact-match seeds. Thus, it would be interesting to compare n = 1 to n = 2 using optimized spaced/subset seed patterns: this was done previously, and n = 2 was worse [START_REF] Frith | Improved search heuristics find 20 000 new alignments between human and mouse genomes[END_REF].

Sparsity via words

Let us now see how seeds starting with a compare to every-4th seeding. For a given seed length, the random hit counts are the same (as expected), but seeds starting with a have lower sensitivity (Figure 2A). This is not too surprising, because every-4th seeding is sparse in just one sequence, but seeds starting with a are sparse in both sequences. Seeds starting with ry also have the same random hit counts, and their sensitivity is closer to (but still less than) that of every-4th seeds. On the other hand, seeds starting with rr have worse sensitivity. This supports the idea that non-overlapping words are good and highly-overlapping words are bad. A number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q a a a a a a

B

number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q every-8th words: abb words: avv

C

number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q every-10th words: abbb words: avvv number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q every-13th words: abbbb words: avvvv 5 Seeds starting with abb have a sparsity factor of 4 3 /3 2 ≈ 7.1, and they perform slightly worse than every-8th seeds (Figure 2B). On the other hand, they perform better than seeds starting with avv (v = any base except t). Seeds starting with abbb (sparsity 9.5) or abbbb (sparsity 12.6) show similar results (Figure 2C-D), confirming the advantage of non-overlapping words.

Minimal-variance words

We can perhaps do better by using longer words with some overlap. Seeds starting with ry are the same as seeds starting with ryn (where n is any base), so it may be better to replace ryn with a less-overlapping set of length-3 words.

It is not obvious how best to quantify "amount of overlap", but one idea is to use variance of occurrence number in random sequences. Let us try these two measures of overlap: VMR1 (variance-tomean ratio from Equations 7,8) and VMR2 (from Equations 9,10).

It is also unclear how to find a set of words that minimizes VMR1 or VMR2, because the number possible sets is enormous. Brute-force search is feasible if we restrict ourselves to a 2-letter ry alphabet.

Such words can indeed boost sensitivity. For example, the words rrry,ryrr,ryyr,yyyr have lower VMR2 than rynn (Table 3), and seeds starting at these words have better sensitivity (Figure 3A). We can do better still with eight length-5 words (Table 3, Figure 3A).

For even longer words, our brute-force search was too slow, so we switched to a heuristic search method (simulated annealing) that does not guarantee to find the minimum possible VMR. This found a set of sixteen length-6 words (Table 3), with even better sensitivity (Figure 3A).

On the other hand, we also found cases where words with lower VMR1 and VMR2 have worse sensitivity (see the Supplement). Thus, a better criterion for choosing a set of words is still to be designed.

Figure 3B-C shows other examples where longer minimum-variance words (from Table 3) improve sensitivity, with no change in sparsity or random hit rate. However, the longer words perform badly when they exceed the seed length (e.g. blue and red points at the right of Figure 3C), showing a danger of too-long words.

Minimizers

We next tested minimizers, with three orderings:

• Alphabetic order.

• cg-order, where cgcg. . . is the minimum sequence. This is representative of methods that have been used in practice [MPB + 17]. • abb-order. This is a novel ordering, inspired by non-overlapping abb. . . words. It compares two suffixes using order a<c<g<t at the first position and t=g=c<a at all subsequent positions. Let us first see the sparsity (average distance between seed start coordinates) of these orderings. Alphabetic minimizers have the lowest sparsity (highest density) for a given window length w, and cg A sensitivity relative to every-4th seeding q q q q q q q q q q q sparsity 4 ry rrry,ryrr,ryyr,yyyr eight 5-mers sixteen 6-mers

B

sensitivity relative to every-8th seeding q q q q q q q q q q q sparsity 8 ryy four 5-mers eight 6-mers sixteen 7-mers

C

sensitivity relative to every-16th seeding q q q q q q q q q q 5 6 7 8 9 10 11 12 13 14 q sparsity 16 rryny four 6-mers eight 7-mers number of seed hits between unrelated random sequences Figure 3: Sensitivity (y-axis) and spurious hit count (x-axis) for exact-match seeds with word-based sparsity. Sensitivity was measured on sequence pairs with PAM distance 20. Seed lengths 5-14 were tested, as shown in panel C. In this figure, the sensitivity is shown relative to every-nth sparsity: (% of related sequence pairs found by word-restricted seeds) / (% of related sequence pairs found by every-nth seeds).

A

minimizer window length average distance between seed positions q q q q q q q q q q q q q q q q q q q q q minimizers: alphabetic minimizers: cg minimizers: abb words: abb words: abbb

B

minimizer window length avg distance between seed positions: observed / expected q q q q q q q q q q q q q q q q q q q q q qq qq qqq q qqq q q q qqq q qq q q qq q q qqq q q qq q q qq qq q q q minimizers have higher sparsity (Figure 4), as reported previously [MPB + 17]. Interestingly, abb minimizers have even higher sparsity for w > 10. Now let us see the sensitivity of these minimizers. Taking alphabetic minimizers as an example, if we increase the window size w without changing the seed length, the sensitivity and random hit rate both decrease (Figure 5), as expected. If we then shorten the seeds, the sensitivity and random hit rate increase. Overall, higher w results in lower sensitivity for a given random hit count.

To fairly compare the three kinds of minimizer, we should compare them using different window sizes that achieve the same sparsity. Based on Figure 4A, alphabetic minimizers with w = 7 are comparable to cg minimizers with w = 6, and alphabetic minimizers with w = 20 are comparable to abb minimizers with w = 16. Comparing them thus, cg and abb minimizers are better than alphabetic minimizers (Figure 5). This supports the idea that higher sparsity for a given w improves homology search, which does not seem to have been clearly shown before.

Minimizers versus words

To fairly compare minimizers with words, we should use minimizer window sizes that produce the same sparsity as the words. Figure 6A,B compares words to minimizers with slightly lower sparsity (higher density), giving an unfair advantage to the minimizers. Seeds starting with a perform worse than alphabetic minimizers for PAM distance 20 (Figure 6A), but better for PAM distance 50 (Figure 6B). On the other hand, seeds starting at non-overlapping (ry) or minimum-variance words perform better than alphabetic or minimizers, at both PAM distances. Figure 6C,D compares words to minimizers with about the same sparsity. Seeds at positions of abbb perform slightly better than alphabetic or abb minimizers. In more detail: for a given seed length, the minimizers have worse sensitivity but slightly better specificity.

Next, we compared word-based seeding to the minimizer scheme of the minimap2 software [START_REF] Li | Minimap2: pairwise alignment for nucleotide sequences[END_REF]. This scheme uses exact-match seeds of length 15, with minimizer window w = 10, and an ordering from a particular hash function applied to each 15mer. The expected density is 2/(w + 1) = 0.1818, but we empirically found a slightly higher density, 0.185-0.188, in both random and real sequences. We compared this to twelve length-6 ry words (density 12/2 6 = 0.1875) that minimize VMR2: rrrrry, rryrrr, rryrry, rryyrr, rryyry, ryryrr, ryyyrr, A number of seed hits between unrelated random sequences sensitivity relative to every-4th seeding a a a a a a a a a a q q q q q q q q q q q q q q q q q q q q a q q words: a words: ry words: eight 5-mers minimizers: w=7, alphabetic minimizers: w=6, cg

B

number of seed hits between unrelated random sequences sensitivity relative to every-4th seeding a a a a a a a a a a q q q q q q q q q q q q q q q q q q q q a q q words: a words: ry words: eight 5-mers minimizers: w=7, alphabetic minimizers: w=6, cg

C

number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q words: abbb minimizers: w=20, alphabetic minimizers: w=16, abb

D

number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q words: abbb minimizers: w=20, alphabetic minimizers: w=16, abb ryyyry, ryyyyr, yryyrr, yryyry, yyyyyr.

For this test, random fragments of size 1000 were drawn from human (GRCh38) chromosome 22, then mutated by the PAM process, and the number of conserved 15-mer seeds was counted. At PAM distance 0, minimap has more seeds, i.e. higher density (Figure 7). Nevertheless, at PAM distance ≥ 1, minimap has fewer conserved seeds.

Unification with subset seeds

A straightforward generalization of subset seeds incorporates word-restricted seeding. Recall that subset seeds allow some mismatches (e.g. a↔g and c↔t) at fixed positions. More generally, we can allow any subset of the 16 possible types of match and mismatch.

An example of such a generalized subset seed pattern is: ANNRYrn@y. This specifies seeds of length 9. Positions with A (in this case, the 1st position) allow a:a matches only. Positions with N allow any match. Positions with n allow any match or mismatch. Positions with R allow purine matches only: a:a or g:g. Positions with r allow purine matches or mismatches: a:a, g:g, a:g, g:a. Positions with Y or y likewise allow pyrimidines (c and t). Finally, positions with @ allow any match or transition.

Such a seed pattern has two important properties: sparsity and weight. Sparsity means rarity of compatible positions in a single random sequence. Weight indicates unlikelihood of a chance match to a compatible position: for example, weight 5 means the same unlikelihood as a length-5 exact match.

The hard problem is to design good seed patterns for finding sequences related by a given PAM distance, transition/transversion bias, etc. Fortunately, the seed design software Iedera already allowed this kind of generalized subset seed [START_REF] Kucherov | A unifying framework for seed sensitivity and its application to subset seeds[END_REF]. Here, we used it to design seeds for PAM=20 and κ=3. To constrain the search space in this preliminary study, we only considered patterns based on ry, i.e. having one R or r and one Y or y, and likewise patterns based on ryy. The only other pattern symbols allowed were N, @, and up to 5 ns. Up to 10 transition-tolerant positions other than n were allowed. The resulting patterns are in Table 4. Many other patterns are equally good; we broke ties by preferring ones that start with RY or RYY.

One notable result is that n positions are useful in the ry-based seeds, but not the ryy-based seeds. 10 weight pattern ry-based seeds: sparsity 4 5 RYNN@@ 6 RY@@@@NN 7 RYN@@@@NN 8 RYN@@@@NNN 9 RYN@@@@@@NNN 10 RYN@@@nnNN@@@NN 11 RYN@@@nnNN@@@NNN 12 RYN@@@@NNnn@@@@NNN 13 RYN@@@@NN@nn@@@NNNN 14 RYN@@@@NN@nn@@@@NNNN@ ryy-based seeds: sparsity 8 5 RYY@@@@ 6 RYYN@@@@ 7 RYYN@@@@@@ 8 RYYNN@@@@@@ 9 RYY@@@@@@@@NN 10 RYY@@@@@@@@NNN 11 RYYN@@@@@@@@NNN 12 RYYN@@@@@@@@@@NNN 13 RYYN@@@@@@@@@@NNNN 14 RYYN@@@@@@@@@@NNNNN This is presumably because n positions make overlapping seeds more independent, but sparser seeds have fewer overlaps. As expected, these seed patterns are good for finding sequences that are related by PAM=20 and κ = 3 (Figure 8).

Discussion

When to use sparsity

The aim of seeding methods is to maximize sensitivity while minimizing computational cost (time and memory). Computational cost has two parts: the cost of finding seed matches (c 1 ) and the cost of processing them (c 2 ). Sparsity need not reduce sensitivity, if the seeds are shortened, but it usually increases random seed hits (i.e. c 2 ) for a given sensitivity (Figure 1). A notable exception is exactmatch seeds and every-nth sparsity with small n (e.g. n 2), which does not increase random hits for a given sensitivity (Figure 1). Typically, however, sparsity is beneficial only when long (or rare) seeds do not sufficiently reduce the computational cost.

When to use every-nth sparsity

Every-nth sparsity has better sensitivity per random hits (c 2 ) than either minimizers or word-restricted seeds, see also [START_REF] Almutairy | Comparing fixed sampling with minimizer sampling when using k-mer indexes to find maximal exact matches[END_REF]. So it should be preferred unless its c 1 is significantly worse. It achieves sparsity in just one of the two sequence datasets being compared, which is appropriate for comparing a huge dataset to a moderate-size dataset, e.g. many DNA reads to a moderate-size genome. It might be appropriate for comparing DNA reads to a human genome.

Sparsity in both datasets, with minimizers or word-restricted seeds, is appropriate for "hugeversus-huge" comparisons. A typical example is aligning DNA reads to each other in order to assemble them, which was a major motivation for minimizers [RHH + 04]. Other examples are searching DNA sequences from unknown organisms against a multi-genome database, or checking if DNA data has contamination from other organisms [START_REF] Steinegger | Terminating contamination: large-scale search identifies more than 2,000,000 contaminated entries in Gen-Bank[END_REF].

Words versus minimizers

This study indicates that seeding at minimallyoverlapping words is superior to minimizers. One caveat -bias due to reduced minimizer density at sequence edges -is addressed in the Supplement, and does not change this conclusion. It is important to note, however, that minimizer schemes are still being optimized [MPB + 17, MDK18]. On the other hand, we have barely begun to optimize wordrestricted seeding.

Compared to word-restricted seeds, a minimizer seed match has an extra contextual requirement. A seed match can be destroyed by a mutation inside the seed: this applies equally to both methods. However, minimizers experience an additional effect: a mutation outside the seed can make that seed position no longer a minimizer. This reduces the sensitivity of minimizers, but increases their specificity, which fits our observations.

Our word-restricted seeding has a potential disadvantage: there is no upper bound on distance between words. The probability of longer distance decreases exponentially in complex sequence, but 11 A number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q sparsity 4 ry exact-match seeds ry transition seeds Iedera ry-based seeds

B

number of seed hits between unrelated random sequences % of related sequence pairs found q q q q q q q q q q q sparsity 8 ryy exact-match seeds ryy transition seeds Iedera ryy-based seeds not in simple sequence such as polypurine tracts or short-period tandem repeats. Pure simple-sequence similarities are typically not wanted, because their significance is hard to assess and they do not reliably indicate homology.

Further advantages of words

Word-restricted seeding has further advantages over minimizers. Firstly, it can be co-designed with subset seeds. Secondly, it seems likely that word positions can be found faster than minimizer positions. Thirdly, word-restricted seeding is more conducive to efficient indexes. Seed matches are usually found with an index data-structure. There are various kinds of index, but they often include a lookup table for any possible DNA sequence of some length d. This table can be reduced (or d increased) with word-restricted seeding, because only a subset of length-d words are ever considered.

Co-designed seed patterns

The sensitivity benefit of spaced and subset seeds can be enhanced by using, instead of one seed pattern, several co-designed patterns [START_REF] Buhler | Designing seeds for similarity search in genomic DNA[END_REF][START_REF] Sun | Choosing the best heuristic for seeded alignment of DNA sequences[END_REF].

Each pattern tends to find similarities that tend to be missed by the other patterns. This idea could be combined with word-restricted seeding. For example, we could use four different patterns, each starting with one of the minimally-overlapping words RRRY, RYRR, RYYR, YYYR (Table 3). Most interestingly, the best set of words may then not be minimally-overlapping ones, but rather words whose overlaps complement the seed patterns.

Open questions

Our study provides a new motivation for the problem of maximizing the number of non-overlapping words. For our purposes, minimally-overlapping words are especially useful, but we remain unsure how best to quantify overlap. Another challenge is how to search a large number of possible word sets for one with low overlap. More generally, we would like to design a set of word-restricted subset-seed patterns.

A further difficulty is how to optimize wordrestricted seeding when the letter frequencies are unequal. In this case, we cannot simply seek an optimal set of n length-k words, because the sparsity is not constant. It is notable, however, that most natural DNA has near-equal frequencies of r and y.

Going further in the direction of empirical data, it might be useful to optimize word-restricted seeding for a particular sequence set (e.g. a genome). Presumably, it is beneficial to use words that are anti-clumped while tending to avoid repetitive sequence. Minimally-overlapping words avoid some kinds of repeat, e.g. homopolymers. Previously, minimizer ordering was defined by frequency in a particular sequence set [CLJ + 14].

Word-restricted seeding requires fast word-finding. Perhaps some word sets are conducive to fast detection, e.g. the words in the last row of Table 3 share a common prefix.

When we use increasingly long minimum-variance words, with fixed sparsity n, the sensitivity might approach that of every-nth seeding (Figure 2, 3). The seed count of every-nth seeding has zero variance: can the words achieve arbitrarily-low variance? If so, they become arbitrarily close to a universal k-mer hitting set. Perhaps optimized minimizers, minimally-overlapping words, and universal k-mer hitting sets will converge.

Figure 1 :

 1 Figure1: Sensitivity (y-axis) and spurious hit count (x-axis) for exact-match seeds with every-nth sparsity. Sensitivity was measured on sequence pairs with PAM distance 20 (left panel) or 50 (right panel). Seed lengths 5-14 were tested, shown in gray in the left panel.

Figure 2 :

 2 Figure 2: Sensitivity (y-axis) and spurious hit count (x-axis) for exact-match seeds with every-nth or word-based sparsity. Sensitivity was measured on sequence pairs with PAM distance 20. Seed lengths 5-14 were tested, as shown in panel A.

Figure 4 :

 4 Figure 4: Sparsity of minimizers, with three orderings. The diagonal gray line in A, and the horizontal gray line in B, show the expected minimizer sparsity (w + 1)/2.

Figure 5 :

 5 Figure5: Sensitivity (y-axis) and spurious hit count (x-axis) for exact-match seeds at minimizer positions. "w" means window length. Seed lengths 5-14 were tested, shown in gray in the left panel.

Figure 6 :

 6 Figure 6: Sensitivity (y-axis) and spurious hit count (x-axis) for exact-match seeds at either word positions or minimizer positions. Seed lengths 5-14 were tested. Sensitivity was measured on sequence pairs with PAM distance 20 (A, C) or 50 (B, D).

Figure 8 :

 8 Figure 8: Sensitivity (y-axis) and random hit count (x-axis) of seeding methods, for sequences with transition/transversion bias (κ = 3) and PAM distance 20. Seed weights 5-14 were tested. "Transition seeds" allow transition substitutions at all positions.

Table 1 :

 1 Parameters of the T92 DNA model

	PAM κ %g+c %identity transitions per
					transversion
	20	1	50	82	0.5
	50	1	50	64	0.5
	20	3	50	83	1.4

Table 2 :

 2 Non-overlapping DNA words. r = {a, g}; y = {c, t}; b = {c, g, t}.

	word	constructed	maximum
	length words	number	number
	2	ry	4	4
	3	abb	9	9
	4	abbb	27	27
	5	abbbb	81	81
	6	abbbbb	243	251

Table 3 :

 3 Variance-to-mean ratios. Bold values are known to be the minimum possible, for that sparsity and word length.

	words	VMR1 VMR2
	sparsity 4
	ry	0.25	0.5
	ryn	0.25	0.417
	rynn	0.25	0.375
	rrry,ryrr,ryyr,yyyr	0.25	0.344
	rrrry,rryrr,ryryr,ryyrr,
	ryyry,ryyyy,yyyrr,yyyry	0.125	0.25
	rrrrry,rryrry,rryryy,ryrrrr,
	ryrrry,ryryry,ryyrrr,ryyrry,
	ryyryr,ryyryy,ryyyry,ryyyyy,
	yryrry,yyyrrr,yyyrry,yyyyry	0.0938 0.188
	sparsity 8
	ryy	0.375	0.625
	rrrry,yrrry,yrryy,yryyy	0.25	0.45
	rrrrry,yrrrry,yrrryr,yrrryy,
	yrryry,yrryyy,yyryry,yyryyy	0.188	0.372
	rrryrrr,rrryryr,rryrryr,rryyrrr,
	rryyrry,rryyryr,rryyyrr,rryyyyr,
	ryryyrr,ryyyryr,ryyyyyr,ryyyyyy,
	yryyryr,yryyyrr,yryyyyr,yyryyrr	0.176	0.329
	sparsity 16
	rryny	0.438	0.688
	rryryy,rryyry,rryyyr,rryyyy	0.375	0.635
	rryryry,rryryyr,rryryyy,rryyryr,
	rryyryy,rryyyry,rryyyyr,rryyyyy	0.281	0.576

Table 4 :

 4 Seed patterns designed by Iedera for PAM 20, κ = 3, alignment length 64.

Acknowledgments

We are grateful to Paul for suggesting seeds starting with a, and Shotaro Tadachi for investigating word-free tracts in human DNA. GK was partially funded by RFBR, project 20-07-00652, and joint RFBR and JSPS project 20-51-50007. LN was partially funded by ANR, ASTER project ANR-16-CE23-0001.


Karel Brinda
email: karel.brinda@hms.harvard.edu

Michael Baym

Gregory Kucherov

Karel Břinda

Simplitigs as an efficient and scalable representation of de Bruijn graphs

Keywords: Sequence analysis, de Bruijn graphs, -mers, representation, scaling k

come

Simplitigs as an efficient and scalable representation of de Bruijn graphs

Background

Advances in DNA sequencing started the golden age of biology when previously unobservable phenomena can be studied on an unprecedented scale. However, sequencing capacity has been growing faster than computer performance and memory, and also faster than available human resources. Today large amounts of sequence data are available. Consequently, traditional sequence-based representations and sequence alignment-based techniques [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF][START_REF] Smith | Identification of common molecular subsequences[END_REF][START_REF] Gotoh | An improved algorithm for matching biological sequences[END_REF] have become less suitable for real-life scenarios due to the space-and time-complexities they impose and their inefficiency in handling polymorphisms.

De Bruijn graphs provide an elegant solution for genomic data representation. They build on top of the concept of k -mers, which are substrings of a fixed length of the genomic strings to be represented, such as sequencing reads, k genomes, or transcriptomes. For a given -mer set, the corresponding de Bruijn graph is a directed graph with the k -mers being vertices and long overlaps between pairs of these -mers indicating edges (Methods). There is k k 1 k an obvious correspondence between -mer sets and de Bruijn graphs, and we can use both terms interchangeably. k

If is chosen appropriately, de Bruijn graphs capture substantial information about the sequenced molecules as k these correspond to some walks in the graph.

The use of de Bruijn graphs is ubiquitous in sequence analysis. Genome assembly uses the property that sequenced molecules form paths [START_REF] Idury | A New Algorithm for DNA Sequence Assembly[END_REF][START_REF] Pevzner | 1-Tuple DNA Sequencing: Computer Analysis[END_REF][START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF] , which is exploited in numerous modern assemblers [START_REF] Zerbino | Velvet: algorithms for de novo short read assembly using de Bruijn graphs[END_REF][START_REF] Simpson | ABySS: A parallel assembler for short read sequence data[END_REF][START_REF] Bankevich | SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing[END_REF][START_REF] Chikhi | Space-efficient and exact de Bruijn graph representation based on a Bloom filter[END_REF][START_REF] Li | MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF][START_REF] Souvorov | SKESA: strategic k-mer extension for scrupulous assemblies[END_REF] . On the other hand, alignment-free sequence comparison follows the idea that similar sequences share common -mers, and k comparing de Bruijn graphs provides thus a good measure of sequence similarity [START_REF] Zielezinski | Alignment-free sequence comparison: benefits, applications, and tools[END_REF] . This involves applications of de Bruijn graphs to variant calling and genotyping [START_REF] Iqbal | De novo assembly and genotyping of variants using colored de Bruijn graphs[END_REF][START_REF] Bradley | Rapid antibiotic-resistance predictions from genome sequence data for Staphylococcus aureus and Mycobacterium tuberculosis[END_REF][START_REF] Shajii | Fast genotyping of known SNPs through approximate k-mer matching[END_REF][START_REF] Sun | Toward fast and accurate SNP genotyping from whole genome sequencing data for bedside diagnostics[END_REF][START_REF] Nordström | Mutation identification by direct comparison of whole-genome sequencing data from mutant and wild-type individuals using k-mers[END_REF] , transcript abundance estimation [START_REF] Bray | Near-optimal probabilistic RNA-seq quantification[END_REF] , and metagenomic classification [START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF][START_REF] Ames | Scalable metagenomic taxonomy classification using a reference genome database[END_REF][START_REF] Břinda | ProPhyle: An accurate, resource-frugal and deterministic DNA sequence classifier[END_REF][START_REF] Corvelo | taxMaps: comprehensive and highly accurate taxonomic classification of short-read data in reasonable time[END_REF] . In the latter application, -mer-based classifiers perform best among all classifiers in k inferring abundance profiles [START_REF] Ye | Benchmarking Metagenomics Tools for Taxonomic Classification[END_REF] , which suggests that de Bruijn graphs truthfully approximate the graph structures of bacterial pan-genomes, even if constructed from noisy assemblies from incomplete databases. Even if more advanced pan-genome graph representations are available, such as variation graphs [START_REF] Paten | Genome graphs and the evolution of genome inference[END_REF] , de Bruijn graphs with large -mer lengths are still useful for indexing [START_REF] Sirén | Indexing Variation Graphs[END_REF][START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF] . k 2 As de Bruijn graphs are one of the primary data structures in much of sequence analysis, the efficiency of many algorithms is directly tied to the efficiency of computation and representation of the graph. De Bruijn graphs can be readily computed through a scan of the datasets including the raw reads, genomes or multiple sequence files. In practice, such a scan often consists in -mer counting as this allows efficient denoising of the graph, e.g., by k removing low-frequency -mers corresponding to sequencing errors in the reads. Algorithms for -mer counting k k

have been extensively studied and many well-engineered software solutions are available [START_REF] Melsted | Efficient counting of k-mers in DNA sequences using a bloom filter[END_REF][START_REF] Marçais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF][START_REF] Deorowicz | Disk-based k-mer counting on a PC[END_REF][START_REF] Rizk | DSK: k-mer counting with very low memory usage[END_REF][START_REF] Roy | Turtle: Identifying frequent k-mers with cache-efficient algorithms[END_REF][START_REF] Crusoe | The khmer software package: enabling efficient nucleotide sequence analysis[END_REF][START_REF] Deorowicz | KMC 2: fast and resource-frugal k-mer counting[END_REF][START_REF] Kokot | KMC 3: counting and manipulating k-mer statistics[END_REF][START_REF] Pandey | Squeakr: an exact and approximate k-mer counting system[END_REF] .

On the other hand, de Bruijn graph representations have received much less attention. The most commonly used representation are unitigs, which are strings resulting from compaction of -mers along maximal paths with k non-branching nodes [START_REF] Chikhi | On the Representation of De Bruijn Graphs[END_REF][START_REF] Chikhi | Compacting de Bruijn graphs from sequencing data quickly and in low memory[END_REF] . Unitigs have many advantages: the representation is "textual", in the form of a set of sequences that contain each -mer exactly once while preserving graph topology. However, unitigs impose large k resource overhead for many types of de Bruijn graphs and do not scale well when a lot of variation is included.

Specifically, with a high proportion of branching nodes, unitigs become fragmented, in an extreme case up to the level of individual -mers. Subsequently, unitig computation and storage may require inappropriately large k resources and become prohibitive in variation-heavy applications, e.g., in bacterial pan-genomics.

In this paper, we propose simplitigs as a compact, efficient and scalable representation of de Bruijn graphs.

Simplitigs generalize the unitig representation by relaxing the restriction of stopping at branching nodes. We present an algorithm for rapid simplitig computation from a -mer set and implement it in a tool called ProphAsm. k

ProphAsm proceeds by loading a -mer set into memory and a greedy enumeration of maximal vertex-disjoint k paths in the associated de Bruijn graph. We use ProphAsm to evaluate the improvement of simplitigs over unitigs, in terms of two main characteristics: the cumulative sequence length (CL) and the number of sequences (NS). We demonstrate that greedily computed simplitigs are close to theoretical bounds in practical applications and provide, compared to unitigs, a substantial improvement in memory requirements and speed in applications such as -mer k matching.

3

Results

The concept of simplitigs

We developed the concept of simplitigs to efficiently represent de Bruijn graphs of sequence data ( Fig. 1 ).

Simplitigs are a generalization of unitigs and correspond to spellings of vertex-disjoint paths covering a given de Bruijn graph ( Fig. 1a , Methods). Consequently, maximal simplitigs are such simplitigs where no two simplitigs can be merged by a overlap (Methods). Note that unitigs and -mers are also simplitigs, but not maximal in k )

( 1 k general ( Fig. 1b ). The main conceptual difference between maximal simplitigs and maximal unitigs is that simplitigs are not limited by branching nodes, which allows for further compaction, with a benefit increasing proportionally to the amount of branching nodes in the graph.

To compare simplitig and unitig representations, we created a benchmarking procedure based on the two characteristics: the number of sequences (NS) and their cumulative length (CL) (Methods, example in Fig. 2 ).

While NS determines the number of records to be kept in memory, CL largely determines the total memory needed.

NS and CL are readily bounded from below by one and by the number of -mers, respectively, and they are also k tightly connected ((eq 1) in Methods). As every step of compaction decreases both NS and CL ( Fig. 1b , Methods), we can optimize them jointly. However, finding an optimal simplitig representation translates to the vertex disjoint path coverage problem. While NP-hard for general graphs (by reduction from the well-known NP-hard problem of computing a Hamiltonian path), the problem may be tractable for observed de Bruijn graphs (Methods).

Since practical applications do not require optimal simplitigs, we prioritized speed and designed a greedy algorithm for their rapid computation ( Alg. 1 , Methods). In an iterative fashion, the algorithm selects an arbitrary -mer as a k seed of a new simplitig and keeps extending it forwards and then backwards as long as possible, while removing the already used -mers from the set. This process is repeated until all -mers are covered. Loading -mers into k k k memory and simplitig computation are linear in the length of the input and the number of -mers, respectively, k 4 and the memory footprint is linear in the number of -mers. We implemented Alg. 1 in a program called k

ProphAsm, available at https://github.com/prophyle/prophasm .

Simplitigs of model organisms

We first evaluated simplitig and unitig representations on assemblies of six model organisms ( Fig. 2 ). As different applications of de Bruijn graphs call for different -mer lengths, we sought to characterize the NS and CL scaling k for both representations with growing, as well as the effect of the species' genome size. Therefore, selected model k organisms were evaluated in an increasing order of the genome size and benchmarked for both representations on a range of -mer lengths corresponding to common alignment-free-based applications [START_REF] Bray | Near-optimal probabilistic RNA-seq quantification[END_REF][START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF][START_REF] Břinda | Rapid inference of antibiotic resistance and susceptibility by genomic neighbour typing[END_REF] . k

We observe that simplitigs always provide substantially better performance than unitigs ( Fig. 2 ). In particular, they quickly approach the theoretical lower bounds for both characteristics tested. Every data set has a range of -mer k lengths where the difference between simplitigs and unitigs is very large, and after a certain threshold, the difference almost vanishes. While for short genomes this threshold is located at smaller -mer lengths than those k typically used in alignment-free applications (e.g., for E. coli ), for bigger genomes this threshold has not 7 k ≈ 1 been attained on the tested range and seems to be substantially shifted towards large -mers (e.g., B. mori ). k

Interestingly, maxima of the NS and CL values for both representations occur very close to the value , G k = log 4 where is the genome size ( Fig. 2 ). This is readily explained by edge saturation: for values of up to , an

G k G log 4
overwhelming fraction of all -mers belong to the genome, which makes the de Bruijn graph branch at nearly 4 k k every node. As a consequence, unitigs are then essentially reduced to individual -mers and their number grows k exponentially whereas simplitigs stay compact on the whole range of -mer lengths. Starting from the k G k = log 4 graph starts to form longer non-branching paths, which drives down the NS and CL of unitigs, and they approach those of simplitigs. However, the difference between simplitigs and unitigs in their count and length may stay considerable even for larger values of , especially in case of large eukaryotic genomes. 8

Performance assessment

We then analyzed the computational resources that had been used to compute the simplitigs using ProphAsm and the unitigs using BCALM ( Fig. 3 ). Both representations were computed in the environment of a computational cluster, with individual experiments deployed as parallel jobs using SLURM (Methods). Even though we primarily focused on the total CPU time, we also tested BCALM using 4 threads to evaluate the effect of parallelization.

We observe that throughout our experiments ProphAsm outcompeted BCALM, both in terms of the memory consumption and the CPU time, even if BCALM was run with 4 threads ( Fig. 3 ). The only exception was the memory consumption for H. sapiens , suggesting that BCALM may be more memory-efficient for large genomes.

Importantly, ProphAsm used a consistent amount of resources: memory 38-51 bytes of RAM per -mer in the k dataset, a limited CPU time and no additional disk space. On the other hand, BCALM's resource usage was less predictable and less consistent across experiments and required frequent trial-and-error resource adjustments and re-running (Methods).

Besides the comparatively high memory and CPU time requirements, the most challenging resource was the disk space. In order to fit within the available disk capacity, we reduced the BCALM consumption using the '-maxdisk' parameter. However, despite that we requested BCALM to use no more than 30 GB per experiment, a random manual inspection revealed a substantially higher use -for instance, 116 GB of disk space for H. sapiens with and 4 threads. The extensive disk space consumption also disabled a similar comparison on a desktop

7 k = 1 computer.
Overall, the better performance of ProphAsm can be explained by simplitigs being fundamentally easier to compute than unitigs and by BCALM being optimized for a particular class of de Bruijn graphs that emerge in assembly-like applications. As the ProphAsm resource usage is highly predictable, it appears to be more suitable for the use on clusters in many parallel instances, unless the number of -mers in the dataset exceeds a critical threshold k determined by the available RAM.

Simplitigs of bacterial pan-genomes

We then sought to evaluate the impact of additional variation in a de Bruijn graph ( Fig. 4 ). Such variation may originate in polymorphisms, varying gene content in a population of genomes that are represented jointly, in haplotypes of viral quasispecies, or in sequencing errors in case of graphs constructed directly from sequencing reads. In all these cases, many nodes of the de Bruijn graph become branching and new paths emerge. To model gradually increasing variation, we used bacterial pan-genomes with different levels of sampling. Given the high diversity and variability of bacteria, de Bruijn graphs provide a convenient option for computational pan-genomes [START_REF] Marschall | Computational pan-genomics: status, promises and challenges[END_REF] . Such pan-genomes can be constructed from draft assemblies or even directly from sequencing reads, and thanks to bacterial genomes being short and haploid, the information captured by the graphs is sufficient for many analyses.

We first constructed a pan-genome of N. gonorrhoeae , and characterized unitigs and simplitigs as a function of pan-genome size ( Fig. 4a ). We used 1,102 draft assemblies of clinical isolates from the Prevention's Gonococcal Isolate Surveillance Project [START_REF] Grad | Genomic Epidemiology of Gonococcal Resistance to Extended-Spectrum Cephalosporins, Macrolides, and Fluoroquinolones in the United States, 2000-2013[END_REF] , from which we built a series of de Bruijn graphs using an increasing number of genomes. Consistent with previous experiments ( Fig 2ab , ), both representations perform comparably well

1 k = 3
when only one bacterial genome is included ( Fig. 4a ). However, as the number of genomes or -mers grows, the k NS and CL grow as well, but with an increasing gap between unitigs and simplitigs; importantly, the latter stay close to the theoretical lower bounds. When the pan-genome size is measured via the number of genomes included, the CL and NS resemble logarithmic functions for both unitigs and simplitigs ( Fig. 4a , left-hand column). However, when the number of -mers included is used instead, the NS and CL functions act as affine functions ( Fig. 4a , k right-hand column). This suggests that a pan-genome -mer count and a species-specific slope may be used as the k predictors of simplitig performance in future applications.

To analyze the relative benefit of simplitigs with growing de Bruijn graphs, we evaluated the NS and CL reduction ratio of simplitigs over unitigs in different configurations ( Fig. 4b ). We used the same N. gonorrhoeae dataset and considered also another dataset of S. pneumoniae , consisting of 616 draft Illumina assemblies of isolates from a carriage study of children in Massachusetts, USA [START_REF] Croucher | Population genomics of post-vaccine changes in pneumococcal epidemiology[END_REF][START_REF] Croucher | Population genomic datasets describing the post-vaccine evolutionary epidemiology of Streptococcus pneumoniae[END_REF] . For both species and for , we constructed a 8, 1 k = 1 3 series of de Bruijn graphs as previously, but this time we visualized the NS and CL reduction ratios. In all cases, the NS reduction ratios eventually stabilized at values close to 3, following an L-shape ( ) or being almost constant 8 k = 1

(
). The CL reduction ratio admitted approximately a logarithmic dependence on the number of genomes and

1 k = 3
still resembled a linear dependence on the number of -mers. Overall, these experiments provided further evidence k that the benefit of simplitigs over unitigs grows with the increased proportion of branching nodes in a de Bruijn graph or with increasing data in case of pan-genome reference structures.

Application of simplitigs for k -mer search

Finally, we sought to demonstrate the benefit of simplitigs in a real application. A major use of de Bruijn graphs consists in -mer matching, which requires the graphs to act as a membership data structure. As both simplitigs k and unitigs are text-based representations, -mer queries can be implemented using an arbitrary full-text k index [START_REF] Mäkinen | Genome-Scale Algorithm Design[END_REF] , notably a Burrows-Wheeler Transform index [START_REF] Ferragina | Opportunistic data structures with applications[END_REF] (sometimes referred to as an FM-index). Here, we used the index of BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] , as one of the best-engineered solutions available, to analyze the impact of unitig replacement by simplitigs.

Single pan-genome

We first evaluated the simplitig improvement on the same N. gonorrhoeae pan-genome ( Fig. 5 ). We considered four different -mer lengths and for each of them, we built three pan-genome representations from k 9, 3, 7, 1 k = 1 2 2 3 the original draft genome assemblies: first, we merged the assemblies as the most straightforward approach to collect all -mers; then we computed pan-genome representations by unitigs and simplitigs. As all the three k representations carry the same -mer set, a full-text index built upon them should provide the same results, but k with a performance reflecting the differences in NS and CL.

We analysed the NS and CL characteristics of the computed representations ( Fig. 5a ). Simplitigs improved NS and CL over unitigs by factors of 3.1×-3.2× and 1.5×-1.6×, respectively, consistent with Fig. [START_REF] Idury | A New Algorithm for DNA Sequence Assembly[END_REF] . This suggests that memory required for unitigs would be approximately by 50% higher than for unitigs across all the values of k considered. We also compared simplitigs and unitigs to assemblies; both improved the CL by two orders of magnitude while the NS stayed comparable for simplitigs and increased twofold for unitigs.

We then evaluated -mer query performance and memory footprint ( Fig. 5b ). For each representation, we k constructed a standard BWA index and matched 10 million random -mers from the pan-genome using BWA k fastmap [START_REF] Li | Exploring single-sample SNP and INDEL calling with whole-genome de novo assembly[END_REF] , with no restriction on maximum size of the suffix-array interval to ensure evaluation correctness (Methods). Surprisingly, simplitigs improved memory by factors of 2.7× -3.0× ( Fig. 5b ), thus twice as much as we 13 previously anticipated. This is explained by the fact that the underlying full-text engine has to keep information about individual sequences in memory as separate records. As NS grows, it has a negative impact on both the memory footprint and query speed. Nevertheless, since simplitigs provided a substantial reduction in NS over unitigs, this overhead has been reduced. Therefore, the excessive number of unitigs observed throughout our experiments ( Fig. 1 and Fig. 2 ) provides a further argument for replacing unitigs by simplitigs when possible. 14

Multiple pan-genomes

We evaluated the performance of the simplitig representation for simultaneous indexing of a large number bacterial pan-genomes ( Fig. 6 ). We downloaded all complete bacterial genomes from Genbank that had not been excluded

from RefSeq (as of May 2020; 9,869 records out of which 9,032 had genomic sequences available; Methods). We restricted ourselves to complete genomes as draft genomes in Genbank are substantially impacted by false genetic variability [START_REF] Merchant | Unexpected cross-species contamination in genome sequencing projects[END_REF][START_REF] Lu | Removing contaminants from databases of draft genomes[END_REF][START_REF] Steinegger | Terminating contamination: large-scale search identifies more than 2,000,000 contaminated entries in GenBank[END_REF] that is particularly common in bacterial studies, mainly due to the contaminant DNA [START_REF] Goig | Contaminant DNA in bacterial sequencing experiments is a major source of false genetic variability[END_REF] . By grouping individual genomes per species, we obtained 3,179 bacterial pan-genomes which we call the "All" dataset.

After computing simplitigs and unitigs per species, we merged the obtained representations and constructed indexes using BWA; all this was done for to evaluate the impact of the -mer length. As the unitig 9, 3, 7, 1 k = 1 2 2 3 k index could not fit into RAM of our desktop computer for any , we also created the "Solid" dataset by omitting k pan-genomes with less than 11 genomes; this resulted in 112 pan-genomes with 3,958 genomes. We provide all the constructed pan-genomes in the form of simplitigs on Zenodo ( 10.5281/zenodo.3800713 ).

First we analyzed the obtained simplitig and unitig representations of both datasets ( Fig. 6a ). We observe that simplitigs provide substantial improvement in both test characteristics. In the Solid dataset, NS and CL were reduced by simplitigs by factors of 3.1-4.5 and 1.4-1.9, respectively; and in the All dataset, NS and CL were reduced by factors of 3.0-4.3 and 1.2-1.4, respectively; all consistent with the scaling observed previously ( Fig. 4 , Fig. 5 ).

While the improvement in NS was almost identical in both datasets (consistent with the top-right graph in Fig 4b ), the improvement in CL was clearly better in the Solid dataset. Indeed, as the vast majority of pan-genomes in the All dataset contained only one genome, the de Bruijn graphs had a comparatively low number of branching nodes, therefore the difference between simplitigs and unitigs was less striking (consistent with the values for small pan-genome sizes in Fig 4b ). We also observe that, in contrast to unitigs, -mer length had only little impact on k the CL of simplitigs within the tested range, which provides better guarantees on required computational resources in future applications.

We then measured the performance of -mer lookup ( Fig. 6b ). Both on a desktop and on a cluster, we evaluated k memory footprints, index loading time and time to match ten million random -mers from the index using BWA k 15 (Methods). We observed that simplitigs substantially improved the memory footprint and index loading times. For , simplitigs largely improved the matching times, where the difference is caused by "ghost -mers" on unitig 9 k = 1 k borders; these were more common in this experiment due to the short -mer length and the high number of k unitigs. For higher -mer lengths, simplitigs still provided a moderate improvement in the matching rate. We note k that the query time with BWT-based -mer indexes is dominated by high-frequency -mers. As these are equally k k frequent in simplitigs and unitigs, the observed performance is similar unless many ghost -mers emerge on k sequence borders as seen previously. On the desktop, the unitig All-index could not be evaluated as it did not fit into memory, and the outlier for may be the result of memory swapping. 9 k = 1 17

Discussion

We introduced the concept of simplitigs, a generalization of unitigs, and demonstrated that simplitigs constitute a compact, efficient and scalable representation of de Bruijn graphs for various types of genomic datasets. The two representations share many similarities: they are text-based and individual strings correspond to spellings of vertex-disjoint paths. Both representations can be seen as irreversible transforms, taking a set of input strings and producing a new set of strings preserving the -mer sets. In both cases, the resulting files can easily be manipulated k using common Unix tools, compressed using standard compression techniques, and indexed using full-text indexes.

The main difference consists in that simplitigs do not explicitly carry information about the topology of the de Bruijn graph. Also simplitigs are not expected to have direct biological significance -neighboring segments of the same simplitig may correspond to distant parts of the same nucleic acid or even to different ones. Nevertheless, unitigs can always be recomputed from simplitigs, but this step is not required for many common applications.

Furthermore, a concept analogous to simplitigs, called disjointigs, was recently introduced in the context of genome assembly using A-Bruijn graphs [START_REF] Lin | Assembly of long error-prone reads using de Bruijn graphs[END_REF][START_REF] Kolmogorov | Assembly of long, error-prone reads using repeat graphs[END_REF] , suggesting that simplitigs may be useful beyond the context of topology-oblivious applications.

We provided ProphAsm, a tool implementing a greedy heuristic to efficiently compute maximal simplitigs from a k -mer set. ProphAsm is a spin-off of the ProPhyle software ( https://prophyle.github.io , [START_REF] Břinda | ProPhyle: An accurate, resource-frugal and deterministic DNA sequence classifier[END_REF][START_REF] Břinda | Novel computational techniques for mapping and classification of Next-Generation Sequencing data[END_REF] ) for metagenomic classification, allowing efficient indexing of -mers propagated to individual nodes of a phylogenetic tree. k ProphAsm presents a "naive" implementation of the greedy heuristic ( Alg 1) that can be further improved. For instance, a hash table with better memory management may reduce the memory requirement by a factor of 2.5 [START_REF] Li | Revisiting hash table performance[END_REF] .

Additional memory reduction may be also achieved in a similar fashion as previously done for unitigs [START_REF] Chikhi | Compacting de Bruijn graphs from sequencing data quickly and in low memory[END_REF][START_REF] Pan | Fast de Bruijn Graph Compaction in Distributed Memory Environments[END_REF][START_REF] Guo | deGSM: memory scalable construction of large scale de Bruijn Graph[END_REF] .

Nevertheless, on the studied data, ProphAsm strikingly outcompeted BCALM in all characteristics measured, with the only exception of memory in the case of H. sapiens .

The observed lower requirements of ProphAsm can be attributed to two major differences between unitigs and simplitigs. First, unitigs, and BCALM as one of the reference programs, are designed for assembly applications, where only a small proportion of nodes is usually branching. Moreover, experiments are usually run on a cluster in series, with multiple threads per experiment, rather than in parallel, and the tools can use extensive disk space.

However, when unitigs are tested as a general representation with restricted disk space, the required resources grow substantially. Second, simplitigs are likely to be, from a computational perspective, fundamentally easier than unitigs.

A challenging but also promising feature of simplitig representation is the ambiguity of maximal simplitigs. This is in sharp contrast to maximal unitigs, which are uniquely defined (up to the order, reverse complementing, and cycles). In practice, every algorithm for simplitig computation has to decide which edge will be included at each branching node. Here, we prioritized speed, simplitigs were constructed progressively and lexicographically minimal edges were used in a case of ambiguity. Therefore, the final maximal simplitigs were only dependent on the choice of seeding -mers; these are determined by the specific implementation of 'std::unordered_set' in the C++ k standard library. Nevertheless, characteristics other than speed could readily be prioritized instead. For instance, a more sophisticated heuristic may shift the CL and NS closer to the optimum. One could also aim at adding biological meaning to simplitigs, e.g., by preferring those paths that are better supported by sequence data. On the other hand, streaming algorithms for operations such as merging or intersecting may require specific prescribed forms of unitigs. Finally, simplitigs may be optimized for entropy in order to maximize their compressibility.

The data presented in this paper highlight the scaling of computational resources as more sequencing data become available. The studied N. gonorrhoeae dataset constitutes a relatively complete image of a bacterial population in a geographical region at a given time scale. As such, it can be used to model the "state of completion" of -mer-based k pan-genome representations. On the other hand, the multiple pan-genomes experiment shows how simplitigs perform when a large number of such pan-genomes, although in different states of completion, are considered simultaneously and queried using a single BWT index. Overall, the presented experiments allow us to predict the resources for species for which only limited sequence data are available at present, but more are likely to be generated in the future. Importantly, with more data available, comparative benefits of simplitigs over unitigs grow as we have shown throughout the paper. As the growth of public databases negatively impacts the accuracy of those algorithms that operate on top of lossily represented de Bruijn graphs [START_REF] Nasko | RefSeq database growth influences the accuracy of k-mer-based lowest common ancestor species identification[END_REF] , simplitigs provide a promising solution offering both exactness and scalability.

In modern bioinformatics applications multiple de Bruijn graphs are often considered simultaneously; the resulting structure is usually referred to as a colored de Bruijn graph [START_REF] Iqbal | De novo assembly and genotyping of variants using colored de Bruijn graphs[END_REF] and the associated data structures have been widely

studied [START_REF] Bowe | Succinct de Bruijn Graphs[END_REF][START_REF] Holley | Bloom Filter Trie: an alignment-free and reference-free data structure for pan-genome storage[END_REF][START_REF] Solomon | Fast search of thousands of short-read sequencing experiments[END_REF][START_REF] Muggli | Succinct colored de Bruijn graphs[END_REF][START_REF] Sun | AllSome Sequence Bloom Trees[END_REF][START_REF] Pandey | Mantis: A Fast, Small, and Exact Large-Scale Sequence-Search Index[END_REF][START_REF] Yu | SeqOthello: querying RNA-seq experiments at scale[END_REF][START_REF] Almodaresi | A space and time-efficient index for the compacted colored de Bruijn graph[END_REF][START_REF] Harris | Improved representation of sequence Bloom trees[END_REF][START_REF] Holley | Bifrost -Highly parallel construction and indexing of colored and compacted de Bruijn graphs[END_REF][START_REF] Bradley | Ultrafast search of all deposited bacterial and viral genomic data[END_REF][START_REF] Bingmann | COBS: A Compact Bit-Sliced Signature Index[END_REF] . Although we touched upon this setting in the Multiple pan-genomes section, exploiting the similarity between individual de Bruijn graphs for further compression in simplitig-based approaches is to be addressed in future work. In many applications, including some of the traditional alignment-free methods [START_REF] Zielezinski | Alignment-free sequence comparison: benefits, applications, and tools[END_REF][START_REF] Zielezinski | Benchmarking of alignment-free sequence comparison methods[END_REF] , it is also desirable to consider -mers with counts. In the context of de Bruijn graphs, this leads to the so-called k weighted variant of the problem [START_REF] Pandey | deBGR: an efficient and near-exact representation of the weighted de Bruijn graph[END_REF] . The fact that frequencies of overlapping -mers are usually similar suggests k that -mers can be grouped based on frequencies and simplitigs constructed per group . k

Independently and simultaneously with the work we present here, the simplitig representation was recently studied in [START_REF] Rahman | Representation of k-mer sets using spectrum-preserving string sets[END_REF] under the name "spectrum-preserving string sets". The associated UST tool follows a similar greedy strategy to ProphAsm, although operating on unitigs constructed by BCALM rather than on the original -mers. As we k demonstrated throughout this paper, unitigs are prohibitive for highly branching de Bruijn graphs, where a simplitig construction through unitigs may create a burden on resources and easily become intractable. The paper presents a tighter lower bound on the cumulative length of the representation (CL), termed weight, taking into account the graph topology but requiring a computational overhead. The authors also studied the compression properties of simplitigs when combined with standard compression algorithms. On the other hand, the paper does not study the number of sequences (NS) and continuous scaling for parameters such as -mer length or amount of k additional variation included. 20

Conclusions

In this paper, we addressed the question of efficient and scalable representation of de Bruijn graphs. We showed that the state-of-the-art unitig representation may require adequately large computational resources, especially when de Bruijn graphs contain many branching nodes. We introduced simplitigs, which provide a more compact replacement in applications that do not require explicit information on the graph topology, such as alignment-free sequence comparison and -mer indexing. We introduced a heuristic simplitig computation and showed on the k examples of model species that unless the genome is large, even a naive implementation outperforms BCALM, the main state-of-the-art tool for unitigs. We then studied applications to bacterial pan-genomics and showed that the utility of simplitigs compared to unitigs grows as more data are available. Finally, we demonstrated on the example of full-text -mer indexing that simplitigs can substantially reduce computational resources and allow k computations in situations where unitigs would bring unaffordable costs when exactness should be preserved.

Our work opens many questions and future directions. The presented algorithm for simplitig computation can be improved, parallelized, de-randomized to ensure reproducibility. We anticipate more theoretical advances in the analysis of the minimum vertex-disjoint path cover problem and better mapping to results from other disciplines such as network sciences. We also anticipate improvements in the heuristic approaches that could simplify and parallelize simplitig computation. The nature of the algorithm suggests that simplitigs might be computed online directly from a stream of data such as sequencing reads. We anticipate better implementations and libraries for simplitigs that can be plugged into standard bioinformatics libraries for various programming languages. Another series of questions is related to low-memory transformations of computed simplitigs that would allow precomputing simplitigs on computer clusters, and tailoring to specific applications on standard computers; this includes decreasing , performing set operations on top of simplitig sets and computing maximal unitigs from k simplitigs. A substantial body of work can be anticipated in the direction of text indexing -we showed that simplitigs can be combined with full-text indexes; however, specialized simplitig indexes exploiting simplitig characteristics are yet to be developed. On the theoretical side, many perspectives are open in the direction of orthogonal algorithmic techniques, such as sketching [START_REF] Rowe | When the levee breaks: a practical guide to sketching algorithms for processing the flood of genomic data[END_REF][START_REF] Elworth | To Petabytes and beyond: recent advances in probabilistic and signal processing algorithms and their application to metagenomics[END_REF] , and in relation to other stringology concepts, such as 21 minimal absent words [START_REF] Pinho | On finding minimal absent words[END_REF] and shortest superstrings [START_REF] Gallant | On finding minimal length superstrings[END_REF] . Finally, simplitigs can serve as components in design of various specialized data structures; these can involve membership queries of classes of de Bruijn graphs and colored de Bruijn graphs, where simplitigs may encode not only the -mers themselves, but also additional metadata such k as -mer frequencies or colors. Furthermore, simplitigs will facilitate new full-text-based data structures for k approximate matching, based on the inclusion of the proximity variation, that would be completely intractable with unitigs. Overall, we anticipate that the simplitig representation will become a generic compact representation of de Bruijn graphs, in particular, in the context of large-scale sequence data search engines [START_REF] Bradley | Ultrafast search of all deposited bacterial and viral genomic data[END_REF] and sequence data repositories such as those of NCBI and EBI.

22

Methods

De Bruijn graphs

All strings are assumed to be over the alphabet . A -mer is a string of length . For a string , A, , , } { C G T k k ...s s = s 1 n we define and . For two strings and of length at least , we define the (s)

•• pref k = s 1 • s k (s) •• suf k = s nk+1 • s n s t
(s) (t) suf k = pref k t s→ k k . ⊙ t • suf (t) s k = s |t|k Given a set
of -mers, the de Bruijn graph of is the directed graph with and

K k K V , ) G = ( E V = K
. For every path in , the string is called a spelling of (u, )

| u v}

E = { v ∈ K 2 → k1 v , .., ) p = ( 1 . v p G ⊙ v ⊙ ...⊙ v v 1 k1 2 k1 k1 p
. This definition of de Bruijn graphs is node-centric , as nodes are identified with -mers and edges are implicit. p k Therefore, we can use the terms " -mer set" and "de Bruijn graph" interchangeably. k

Simplitigs

Consider a set of -mers and the corresponding de Bruijn graph . A simplitig graph is a

K k K, ) G = ( E K, ) G′ = ( E′ spanning subgraph of
that is acyclic and the in-degree and out-degree of any node is at most one. It follows from G this definition that a simplitig graph is a vertex-disjoint union of paths, whose spellings we call simplitigs . A simplitig is called maximal if it cannot be extended forward or backward without breaking the definition of simplitig graph. In more detail, a simplitig is maximal if the following conditions hold u ... u

u 1 → k1 2 → k1 → k1 n
• either has no incoming edges in , or for any edge , belongs to another simplitig and it is

u 1 G v, ) ( u 1 ∈ E v not its last vertex,
• either has no outgoing edges in , or for any edge , belongs to another simplitig and it is

u n G u , ) ( n v ∈ E v not its first vertex.
A unitig is a simplitig such that each of the nodes has in-degree

L #kmers (k ) ⋅ N S C = + 1 (eq 1)
As an important consequence, both characteristics are optimized simultaneously.

Greedy computation of simplitigs

The problem of computing maximal simplitigs that are optimal in CL (i.e., also in NS) corresponds to the minimum vertex-disjoint path cover problem [START_REF] Manuel | Revisiting path-type covering and partitioning problems[END_REF] . This is known to be NP-hard in the general case, reducing from the Hamiltonian path problem. However, the complexity for de Bruijn graphs remains an open question. A greedy heuristic to compute maximal simplitigs has been used throughout this paper ( Alg. 1 ). Simplitigs are constructed iteratively, starting from (arbitrary) seeding -mers and being extended greedily forwards and backwards as long k as possible.

ProphAsm implementation

ProphAsm is written in C++ and implements the greedy approach described above ( Alg. 1 ). -mers are encoded K using uint64_t and stored in an std::unordered_map. The choice of extension nucleotides on branching nodes is done based on the lexicographic order. Therefore, the only source of randomness is the choice of seeding -mers by k std::unordered_set::begin; the C++ standard library makes no guarantees on which specific element is considered the first element. ProphAsm does not require any disk space to store intermediate data and its memory consumption corresponded to 38-51 bytes per a unique -mer (in dependence on the allocation), consistent with k [START_REF] Li | Revisiting hash table performance[END_REF] .

Uni-directed and bi-directed models

The uni-directed model, as presented above, is useful for introducing the concepts of unitigs and simplitigs, but is not directly applicable to data obtained using sequencing: since DNA is double-stranded, every string may come from either strand. At the level of -mers, double-strandedness can be accounted for by using canonical -mers, k k

i.e., by pairing-up every -mer with its reverse complement, typically done by taking the lexicographical minimum k of the -mer and its reverse complement. This subsequently requires redefinining de Bruijn graphs to bi-directed k de Bruijn graphs [START_REF] Medvedev | Maximum likelihood genome assembly[END_REF] , which requires a more complex formalism.

Correctness evaluation

The correctness of simplitig computation can be verified using an arbitrary -mer counter. Simplitigs have been k

computed correctly if and only if every -mer is present exactly once and the number of distinct -mers is the same k k as in the original datasets. The correctness of ProphAsm outputs was verified using JellyFish 2 [START_REF] Marçais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF] .

Experimental evaluation -model organisms and performance

Reference sequences for six selected model organisms were downloaded from RefSeq and UCSC Genome Browser:

S. pneumoniae str. ATCC 700669 (accession: NC_011900. Individual experiments were run in parallel on the Harvard Medical School O2 cluster using Snakemake [START_REF] Köster | Snakemake-a scalable bioinformatics workflow engine[END_REF] and SLURM. ProphAsm and BCALM were run with the following parameters, respectively: '-k {kmer-length}' and 25 '-kmer-size {kmer-length} -abundance-min 1 -nb-cores {cores} -max-disk 30000'. As BCALM requires a large undocumented amount of disk space, we used the -max-disk parameter to make a parallel execution of many BCALM jobs feasible. The SLURM specifications of resource allocation for individual species were iteratively adjusted until all jobs would finish; the final required resources are provided in Supplementary Table 1 . Time and memory consumption of jobs were measured independently using GNU Time. Individual jobs were deployed to computational nodes with different hardware configurations, which are specified in Additional File 2 . [START_REF] Croucher | Population genomics of post-vaccine changes in pneumococcal epidemiology[END_REF][START_REF] Croucher | Population genomic datasets describing the post-vaccine evolutionary epidemiology of Streptococcus pneumoniae[END_REF] , and sequenced using Illumina HiSeq)

Supplementary

were downloaded from the SRA FTP server using the accession codes provided in Table 1 in [START_REF] Croucher | Population genomic datasets describing the post-vaccine evolutionary epidemiology of Streptococcus pneumoniae[END_REF] . For each of these datasets, an increasing number of genomes was being taken and merged, and simplitigs and unitigs computed using ProphAsm and BCALM, respectively. This experiment was performed for and To avoid excessive

8 k = 1 1. k = 3
resource usage the functions were evaluated at selected points in an increasing distance: for intervals [10, 100] and [100,+∞] only multiples of 5 and 20 were evaluated, respectively.

Experimental evaluation -full-text k -mer queries

In the single pan-genome experiment, the same 1,102 assemblies of N. gonorrhoeae were merged into a single file.

ProphAsm and BCALM were then used to compute simplitigs and unitigs, respectively, from this file for . Each of the three obtained FASTA files (assemblies, simplitigs, and unitigs) was used to 9, 23, 27, 31 k = 1 construct a BWA index, which was then queried for -mers using 'bwa fastmap -l {kmer-length}'. We used a k modified version of BWA fastmap that reports both the time of index loading and the time of querying ( http://github.com/karel-brinda/bwa , commit e1f907c). Query -mers were generated from the same pan-genome k using WGsim (version 1.10, with the parameters '-h 0 -S 42 -r 0.0 -1 {kmer-length} -N 10000000 -e 0').

For the multiple pan-genome experiment, a list of available bacterial assemblies was downloaded from ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt (2020/05/05). For all assemblies marked as complete (i.e., the "assembly_level" column equal to "Complete genome") and present in RefSeq (i.e., an empty value in the column "excluded_from_refseq"), directory urls and species names were extracted (n=9,869).

These were then used to download the genomes of the isolates using RSync, restricting to genomic sequences only (i.e., files matching '*v?_genomic.fna.gz', n=9,032). The downloaded assemblies were then merged per species in order to collect -mers of individual pan-genomes and used for computing simplitigs and unitigs using ProphAsm k and BCALM, respectively. The obtained simplitig and unitig files were then merged per categories (e.g., simplitigs for k=19) and used to construct a BWA index. The obtained indexes were queried for 1o million -mers using BWA k fastmap as previously. The -mers were generated from the original assemblies of randomly selected 100 genomes k using DWGsim [START_REF] Homer | Whole Genome Simulator for Next-Generation Sequencing[END_REF] (version 0.1.11, with the parameters '-R 0 -e 0 -r 0 -X 0 -y 0 -H -z 42 -m /dev/null -N 10000000 -1 {k} -2 0'); the randomization was performed using 'sort -R'.

Computational setup

The experiments were performed on the HMS O2 research high-performance cluster and on an iMac 4.2 GHz Quad-Core Intel Core i7 with 40 GB RAM. The reproducibility of computation was ensured using BioConda [START_REF] Grüning | Bioconda: sustainable and comprehensive software distribution for the life sciences[END_REF] .

Reads Assemblies

Fig. 1 Alg. 1

 11 Fig. 1 Overview of the simplitig approach. a Textual representations of -mer sets ordered by the degree of k compaction: individual -mers, maximal unitigs, and maximal simplitigs. Every component of a simplitig subgraph k (black arrows) of the de Bruijn graph (all arrows) corresponds to a path, and its spelling constitutes a simplitig (Methods). b Scheme of all possible simplitig representations according to the degree of compaction. While unitigs (dark gray area) correspond to compaction along non-branching nodes in the associated de Bruijn graph, simplitigs (gray area) can also contain branching nodes. Every step of compaction decreases the number of sequences (NS) and their cumulative length (CL) by and by , respectively. Maximal simplitigs may not be determined 1 k 1 uniquely; the simplitig representation can have multiple local optima, depending on which edges were selected at the branching nodes. c The workflow of simplitigs. Simplitigs represent de Bruijn graphs and carry implicitly the same information as unitigs. De Bruijn graphs are usually computed from either assemblies or weighted de Bruijn graphs. Weighted de Bruijn graphs are typically obtained by -mer counting and allow removing noise, e.g., k

k 7 Fig. 2

 72 Fig. 2 Comparison of the simplitig and unitig representations for selected model organisms and a range of -mers. k The number of sequences (NS, millions) and their cumulative length (CL, megabase pairs) for both representations of six model organisms ordered by their genome size: a Streptococcus pneumoniae , 2.22 Mbp; b Escherichia coli , 4.64 Mbp; c Saccharomyces cerevisiae , 12.2 Mbp; d Caenorhabditis elegans , 100 Mbp; e Bombyx mori , 482 Mbp; and f Homo sapiens , 3.21 Gbp. The CL lower bound corresponds to the number of -mers. Full results are available k in Additional File 1 .

Fig. 3

 3 Fig. 3 Comparison of CPU time and memory consumption of ProphAsm and BCALM. Resources to compute unitigs using BCALM (using one and four threads) and simplitigs using ProphAsm (using one thread) of the six model organisms. Full results are available in Additional File 2 .

Fig. 4

 4 Fig. 4 Scaling of simplitigs and unitigs of bacterial pan-genomes as the pan-genome size grows with a better sampling and more within-species variation. a Number of sequences (NS, thousands) and their cumulative length (CL, megabase pairs) for simplitigs (simpl) and unitigs (unit) and the lower bound (lobound) of N. gonorrhoeae and , as a function of the number of genomes (left) and -mers (right, millions) included. b Reduction ratio 1 k = 3 k of simplitigs over unitigs for S. pneumoniae (sp) and N. gonorrhoeae (gc) and , as a function of the 8, 1 k = 1 3 number of genomes (left) and -mers (right, millions) included. Full results are available in Additional File 3 . k

Fig. 5 -ab

 5 Fig. 5 -mer queries for the N. gonorrhoeae pan-genome on top of the draft assemblies, unitigs, and simplitigs. K a Characteristics of the obtained unitigs and simplitigs: number of sequences (NS, thousands) and their cumulative length (CL, megabase pairs). The dot-dash line depicts the CL lower bound corresponding to the number of -mers. k b Time and memory footprint of BWA -mer queries (10 million -mers). Full results including relative k k improvements are available in Additional File 4 .

Fig. 6 -ab

 6 Fig. 6 -mer queries for multiple pan-genomes indexed simultaneously. Bacterial pan-genomes were computed K from the complete Genbank assemblies per individual species. While the All dataset comprises all pan-genomes with no restriction on their size, the Solid dataset comprises only those that contain at least 11 genomes. a Characteristics of the obtained unitigs and simplitigs: number of sequences (NS, millions) and their cumulative length (CL, gigabase pairs). The dot-dash line depicts the lower bound corresponding to the number of -mers. k b Memory footprint, time of index loading and time of matching 10 million -mers using BWA. The bars k correspond to the mean of three measurements (black dots). Full results including relative improvements are available in Additional File 5 .

  nb. of k-mers [M] b nb. of genomes nb. of genomes pan-genome size pan-genome size CL reduction ratio NS reudction ratio nb. of k-mers [M]

Table 1

 1 SLURM resource allocation for ProphAsm and BCALM2 for the performance evaluation.

	Species	ProphAsm (1 core)	BCALM (1 core)	BCALM (4 cores)
		mem	cpu	mem	cpu	mem	cpu
		[GB]	[hours]	[GB]	[hours]	[GB]	[hours]
	S. pneumoniae	10	1	10	1	10	1
	E. coli	10	1	10	1	10	1
	S. cerevisiae	10	1	10	1	10	1
	C. elegans	10	1	10	2	10	2
	B. mori	20	2	30	10	30	2
	H. sapiens	120	4	100	48	100	18
	Experimental evaluation -bacterial pan-genomes

First, 1,102 draft assemblies of N. gonorrhoeae clinical isolates (collected from 2000 to 2013 by the Centers for Disease Control and Prevention's Gonococcal Isolate Surveillance Project

[START_REF] Grad | Genomic Epidemiology of Gonococcal Resistance to Extended-Spectrum Cephalosporins, Macrolides, and Fluoroquinolones in the United States, 2000-2013[END_REF]

, and sequenced using Illumina HiSeq) were downloaded from Zenodo

[START_REF] Grad | Data for "Genomic Epidemiology of Gonococcal Resistance to Extended-Spectrum Cephalosporins, Macrolides, and Fluoroquinolones in the United States, 2000-2013[END_REF]

. Second, 616 draft assemblies of S. pneumoniae isolates (collected from 2001 to 2007 for a carriage study of children in Massachusetts, USA

Acknowledgements

The authors thank Jasmijn Baaijens and Roman Cheplyaka for careful reading and valuable comments, and Kamil Salikhov and Simone Pignotti for helpful discussions at the initial stage of this project. Portions of this research were conducted on the O2 high-performance compute cluster, supported by the Research Computing Group at Harvard Medical School.

Funding KB and MB were partially supported by the David and Lucile Packard Foundation and NIGMS of the National Institutes of Health under award number R35GM133700. GK was partially funded by RFBR, project 20-07-00652, and joint RFBR and JSPS project 20-51-50007.

Availability of data and materials

Additional files

Additional File 1. Detailed information for the single genome experiment: NS, CL and #kmers for unitigs and simplitigs as a function of for the 6 species: a S. pneumoniae , b E. coli , c S. cerevisiae , d C. elegans , e B. mori , k and f H. sapiens .

All benchmarking was performed using ProphAsm 0.1.1 (commit ea28b708) and BCALM 2.2.2 (commit febf79a3). Time and memory footprint were measured using GNU Time.

Declarations Ethics approval and consent to participate

Software

ProphAsm is open source, licensed under the MIT License. The program was developed in C++ and its source code is available from Github ( http://github.com/prophyle/prophasm ). ProphAsm binaries for Linux and OS X are distributed through BioConda [START_REF] Grüning | Bioconda: sustainable and comprehensive software distribution for the life sciences[END_REF] ( https://bioconda.github.io/recipes/prophasm/README.html ). The source code of the version used in this paper was deposited in Zenodo ( 10.5281/zenodo.3887035 ).

Competing interests

No competing interests.

Authors' contributions

KB, MB, GK designed the study, contributed to interpretation of the results, wrote the manuscript, and approved the final manuscript. KB developed the software and performed the data analysis. KB and GK developed the theory.