Compressed Spaced Sux Arrays

Abstract

Spaced seeds are important tools for similarity search in bioinformatics, and using several seeds together often significantly improves their performance. With existing approaches, however, for each seed we keep a separate linear-size data structure, either a hash table or a spaced sux array (SSA). In this paper we show how to compress SSAs relative to normal sux arrays (SAs) and still support fast random access to them. We first prove a theoretical upper bound on the space needed to store an SSA when we already have the SA. We then present experiments indicating that our approach works even better in practice.

-

For the problem of similarity search, we are given two texts and asked to find each suciently long substring of the first text that is within a certain Hamming distance of some substring of the second text. Similarity search has many applications in bioinformatics — e.g., ortholog detection, structure prediction or determining rearrangements — and has been extensively studied (see, e.g., [ 19 ]). Researchers used to first look for short substrings of the first text that occur unchanged in the second text, called seeds, then try to extend these short, exact matches in either direction to obtain longer, approximate matches. This approach is called, naturally enough, “seed and extend”. The substrings’ exact matches are found using either a hash table of the substrings with the right length, or an index structure such as a sux array (SA).

Around the turn of the millenium, Burkhardt and K¨arkk¨ainen [5] and Ma, Tromp and Li [ 15 ] independently proposed looking for short subsequences of the first text that have a certain shape and occur unchanged in the second text, and trying to extend those. A binary string encoding the shape of a subsequence, with 1s indicating positions where the characters must match and 0s indicating positions where they need not, is called a spaced seed. The total number of bits in the binary string is called the seed’s length, and the number of 1s is called its weight. The subsequences’ exact matches are found using either a hash table of the subsequences with the right shape, or a kind of modified SA called a spaced sux array [ 12 ] (SSA).

Burkhardt and K¨arkk¨ainen, Ma et al. and subsequent authors have shown that using spaced seeds significantly improves the performance of seeding and extending. Many papers have been written about how to design spaced seeds to minimize the number of errors (see, e.g., [ 4, 8, 11 ] and references therein), with the specifics depending on the model of sequence similarity and the acceptable numbers of false positives (for which the characters indicated by 1s all match but the substrings are not similar) and false negatives (for which those do not all match but the substrings are still similar) for the application in question. Regardless of the particular application, however, researchers have consistently observed that the best results are obtained using more than one seed at a time. A set of spaced seeds used in combination is called a multiple seed.

Multiple seeds are now a popular and powerful tool for similarity search, but they have a lingering flaw: we keep a hash table or SSA for each seed, and each instance of these data structures takes linear space. For example, SHRiMP2’s [ 7 ] index for the human genome takes 16 GB for each seed. In contrast, Bowtie 2’s [ 14 ] compressed SA for that genome takes only 2.5 GB. This is because a normal SA (which supports only substring matching) can be compressed such that the number of bits per character is only slightly greater than the empirical entropy of the text. Unfortunately, the techniques for compressing normal SAs do not seem to apply directly to SSAs.

In this paper we show how to compress SSAs relative to normal SAs and still support fast random access to them. Whereas the normal SA for a text lists the starting points of the suxes of that text by those suxes’ lexicographic order, the SSA for a text and a spaced seed lists the starting points of the subsequences with the right shape by those subsequences’ lexicographic order. Intuitively, if the seed starts with many 1s, the SSA will be similar to the SA. In Section 2 we formalize this intuition and prove a theoretical upper bound on the space needed to store an SSA when we already have the SA, in terms of the text’s length, the alphabet’s size, and the seed’s length and weight.

In Section 3 we present experiments showing that our approach works even better in practice. That is, even when we implement our data structures using simpler, theoretically sub-optimal components, we achieve better compression than our upper bound predicts. In fact, in practice we can even successfully apply our approach in some cases when the assumptions underlying our upper bounds are violated. However, we still want to improve our compression and random-access times for seeds with low weight-to-length ratios.

We recently learned that Peterlongo et al. [ 17 ] and Crochemore and Tischler [6] independently defined SSAs, under the names “bi-factor arrays” and “gapped sux arrays”, for the special case in which the spaced seed has the form 1a0b1c. Russo and Tischler [ 18 ] showed how to represent such an SSA in asymptotically succinct space such that we can support random access to it in time logarithmic in the length of the text. We note, however, that the spaced seeds used for most applications do not have this form. We also recently learned that Battaglia et al. [ 2 ] used an idea similar to that of spaced seeds in an algorithm for finding motifs with don’t-care symbols. It seems possible our results could be useful in reducing their algorithm’s memory usage. 2

Theory Suppose we want to store an SSA for a text T [0..n 1] over an alphabet of size and a spaced seed S with length ` and weight w. For i < n, let Ti be the subsequence of T [i..n 1] that contains T [j] if and only if i  j and S[j i] = 1. Let Ti0 be the subsequence of T [i..n 1] that contains T [j] if and only if S[j i] = 0. Let SSA be the permutation on {0, . . . , n 1} in which i precedes i0 if either Ti Ti0 , or Ti = Ti0 and T [i..n 1] T [i0..n 1].

For example, if T = abracadabra and S = 101 then

T0 T1 T2 T3 T4 T5 = ar = ba = rc = aa = cd = aa

T6 T7 T8 T9 T10 = db = ar = ba = r = a

T 0

0 T 0

1 T 0

2 T 0

3 T 0

4 T 0 5 = b = r = a = c = a = d

T 0

6 T 0

7 T 0

8 T 0 9 = a = b = r = a and so SSA = [ 10, 3, 5, 7, 0, 8, 1, 4, 6, 9, 2 ], while SA = [ 10, 7, 0, 3, 5, 8, 1, 4, 6, 9, 2 ].

If Ti Ti0 and Ti0 Ti00 , then i precedes i0 in both SSA and SA. In particular, if Ti = Ti0 or Ti0 = Ti00 , then i and i0 have the same relative order in SSA and SA. In our example, T3 = T5 = aa, so 3 precedes 5 in both SSA and SA; T20 = T60 = T90 = a, so 6 precedes 9 and 9 precedes 2 in both SSA and SA.

If we partition SSA into subsequences such that i and i0 are in the same subsequence if and only if Ti = Ti0 , then we can partition SA into the same subsequences. Since there are at most w + w distinct strings Ti, our partitions each consist of at most w + w subsequences. Similarly, if we partition based on Ti0 and Ti00 , then our partitions each consist of at most ` w + ` w subsequences.

For our example, we can partition both SSA and SA into [ 4, 6, 9, 2 ], for Ti0 = a; [ 7, 0 ], for Ti0 = b; [3], for Ti0 = c; [5], for Ti0 = d; [ 8, 1 ], for Ti0 = r; and [ 10 ], for Ti0 = ✏. In this particular case, however, we could just as well partition both SSA and SA into only two common subsequences: e.g., [ 10, 7, 0 ] and [ 3, 5, 8, 1, 4, 6, 9, 2 ].

Consider the permutation SA 1 SSA, which maps elements’ positions in SSA to their positions in SA, and let ⇢ be the minimum number of increasing subsequences into which SA 1 SSA can be partitioned. Since any subsequence common to SSA and SA corresponds to an increasing subsequence in SA 1 SSA, we have ⇢  min( w +w, ` w +` w). In our example, SA 1 SSA = [ 0, 3, 4, 1, 2, 5, 6, 7, 8, 9, 10 ] and ⇢ = 2.

Supowit [ 20 ] gave a simple algorithm that partitions SA 1 SSA into ⇢ increasing subsequences in O(n lg ⇢ ) ✓ O (n min(w, ` w) lg ) time. When applied to SA 1 SSA in our example, Supowit’s algorithm partitions it into [0, 3, 4] and [ 1, 2, 5, 6, 7, 8, 9, 10 ].

Barbay et al. [ 1 ] showed how, given a partition of SA 1 SSA into ⇢ increasing subsequences, we can store it in (2 + o(1))n lg ⇢  (2 + o(1))n min(w, ` w) lg bits and support random access to it in O(lg lg ⇢ ) time. Combining their ideas with later work by Belazzougui and Navarro [3], we can keep the same space bound and improve the time bound to O(1).

To do this, for i  ⇢ , we replace each element in the ith subsequence in SA 1 SSA by a character ai, and store the resulting string R such that we can support random access to it and partial rank queries on it. We then permute R according to SA 1 SSA and store the resulting string R0 such that we can support fast select queries on it. In our example, R = a1a2a2a1a1a2a2a2a2a2a2 and R0 = a1a1a1a2a2a2a2a2a2a2a2.

The partial rank query R.p rank(i) returns the number of copies of R[i] in R[0..i], and the select query R0.selecta(i) returns the position of the ith copy of a in R0. Barbay et al. noted that, for i < n,

(SA 1 SSA)[i] = R0.selectR[i](R.p rank(i)) .

Belazzougui and Navarro showed how we can store R in (1 + o(1))n lg ⇢ bits and support random access to it and partial rank queries on it in O(1) time, and store R0 in (1 + o(1))n lg ⇢ bits and support select queries on it in O(1) time.

In summary, we can store SA 1 SSA in (2 + o(1))n min(w, ` w) lg bits such that we can support random access to it in O(1) time. We will give a longer explanation in the full version of this paper. Since SSA = SA (SA 1 SSA), this gives us the following result: Theorem 2.1 Let T [0..n 1] be a text over an alphabet of size and let S be a spaced seed with length ` and weight w. If we have already stored the sux array SA for T such that we can support random access to SA in time tSA, then we can store a spaced sux array SSA for T and S in (2 + o(1))n min(w, ` w) lg bits such that we can support random access to SSA in tSA + O(1) time. 3

Practice Theorem 2.1 says we can store SSAs for the human Y-chromosome chrY.fa in FASTA format (about 60 million characters over an alphabet of size 5) and SHRiMP2’s three default spaced seeds — i.e., 11110111101111, 1111011100100001111 and 1111000011001101111 — in about 560 MB, in addition to the SA, whereas storing the SSAs na¨ıvely would take about 720 MB. Storing the SSAs packed such that each entry takes dlg 60 000 000e = 26 bits would reduce this to about 580 MB.

To test our approach, we built the SSAs as described in Section 2; computed SA 1 SSA, R and R0 for each SSA; and stored each copy of R or R0 as a wavelet tree. We chose wavelet trees because they are simple to use and often more practical than the theoretically smaller and faster data structures mentioned in Section 2. We ran all our tests described in this section on a computer with a quad-core Intel Xeon CPU with 32 GB of RAM, running Ubuntu 12.04. We used a wavelet-tree implementation from https://github.com/fclaude/libcds and compiled it with GNU g++ version 4.4.3 with optimization flag -O3.

The uncompressed SA took 226 MB, and the six wavelet trees took a total of 215 MB and performed 10 000 random accesses each in 7.67 microseconds per access. That is, we compressed the SSAs into about 60% of the space it would take to store them na¨ıvely and, although our accesses were much slower than direct memory accesses, they were fast compared to disk accesses. Thus, our approach seems likely to be useful when a set of SSAs is slightly larger than the memory and fits only when compressed.

Using the same test setup, we then compressed SSAs for the ten spaced seeds BFAST [9, Table S3] uses for 36-base-pair Illumina reads, which all have weight 18: 1. 111111111111111111 2. 11110100110111101010101111 3. 11111111111111001111 4. 1111011101100101001111111 5. 11110111000101010000010101110111 6. 1011001101011110100110010010111 7. 1110110010100001000101100111001111 8. 1111011111111111111 9. 11011111100010110111101101

Since the first seed consists only of 1s, the SSA we would build for it is the same as the SA. The uncompressed SA again took 226 MB and the 18 wavelet trees for the other nine seeds took a total of 649 MB — so instead of 2.26 GB, we used 875 MB (about 39%) for all ten seeds — and together performed 10 000 random accesses to each of the ten SSAs in about 7 microseconds per access. The left side of the top half of Figure 1 shows how many bits per character (bpc) of the text each SSA took, and the average time per access to each SSA.

We also compressed the SSAs for the ten spaced seeds BFAST uses for 50-basepair Illumina reads, which all have weight 22: 1. 1111111111111111111111 2. 1111101110111010100101011011111 3. 1011110101101001011000011010001111111 4. 10111001101001100100111101010001011111 5. 11111011011101111011111111 6. 111111100101001000101111101110111 7. 11110101110010100010101101010111111 8. 111101101011011001100000101101001011101 9. 1111011010001000110101100101100110100111 10. 1111010010110110101110010110111011 .

Again, the first seed consists only of 1s. This time, the 18 wavelet trees for the other nine seeds took a total of 712 MB; each access took about 8 microseconds. The left side of the bottom half of Figure 1 shows how many bit per character of the text each SSA took, and the average access time per access to each SSA.

If we have a permutation ⇡ 1 on {0, . . . , n 1} stored and ⇡ 2 is any other permutation on {0, . . . , n 1}, then we can store ⇡ 2 relative to ⇡ 1 using the ideas from Section 2. For example, we can store SSAs relative to other SSAs. Suppose we consider the size of each SSA (except the SA) when compressed relative to each other SSA (including the SA), build a minimum spanning tree rooted at the SA, and compress each SSA relative to its parent in the tree. This can reduce our space usage at the cost of increasing the random-access time, as shown for the BFAST seeds on the right side of Figure 1.

There are other circumstances in which we can ignore SSAs’ semantics and consider them only as permutations. For example, spaced seeds can be generalized to subset seeds [ 13 ], such as ternary strings in which 1s indicate positions where the characters must match, 0s indicate positions where they need not, and Ts indicate positions where characters must fall within the same equivalence class (such as the pyrimidines C and T and the purines A and G). It is not dicult to generalize Theorem 2.1 to subset seeds — we will do so in the full version of this paper — but it is also not necessary to obtain practical results. The Iedera tool (available at http://bioinfo.lifl.fr/yass/iedera.php) generates good subset seeds.

A more challenging change is from fixed-length seeds to repetitive seeds [ 12 ]. A repetitive seed is a string in whose repetition the digits indicate which characters must match and how. For example, with respect to the repetitive spaced seed 10110, ATCGATCGGT matches ACCGTTGGGA but not ACCGTTGAGA. Repetitive seeds are useful when looking for approximate matches of substrings that have been extended until they become suciently infrequent. It is not clear how or if we can extend Theorem 2.1 to repetitive seeds. Nevertheless, the LAST tool (available at http://last.cbrc.jp) generates SSAs for repetitive spaced or subset seeds, which we can still try to compress in practice; see also [ 10, 16 ].

Our current goal is to achieve reasonable compression and access times for a set of repetitive subset seeds that we received from Martin Frith, which have average length 19.85 and average weight about 10.44, counting “same equivalence class” digits as 0.5. Unfortunately, at the moment we use nearly 24 bits per entry in the corresponding SSAs (including the overhead for the uncompressed SA), which is only marginally better than the 26 bits we would use with simple packing. Meanwhile, random accesses take about 12 microseconds on average, which is significantly slower than access to a packed array. On the other hand, these seeds have an unusually low average weight-to-length ratio. We used Iedera and LAST to generate SSAs for a set of eight repetitive subset seeds, with average length 17.875 and average weight 12. For these, we used only 20.15 bits per entry, with random accesses taking about 10 microseconds on average. seed Acknowledgments Many thanks to Francisco Claude, Maxime Crochemore, Matei David, Martin Frith, Costas Iliopoulos, Juha K¨arkk¨ainen, Gregory Kucherov, Bin Ma, Ian Munro, Taku Onodera, Gonzalo Navarro, Luis Russo, German Tischler and the anonymous reviewers.

Computer Science, 410:4327–4340, 2009. [3] D. Belazzougui and G. Navarro. Alphabet-independent compressed text indexing. ACM Transactions on Algorithms. To appear. [4] D. G. Brown. Bioinformatics Algorithms: Techniques and Applications, chapter A survey of seeding for sequence alignment, pages 126–152. WileyInterscience, 2008. [5] S. Burkhardt and J. K¨arkk¨ainen. Better filtering with gapped q-grams. Fundamenta Informicae, 56:51–70, 2003. [6] M. Crochemore and G. Tischler. The gapped sux array: A new index structure for fast approximate matching. In Proceedings of the 17th Symposium on String Processing and Information Retrieval (SPIRE), pages 359–364, 2010.

References

[1] J. Barbay , F. Claude , T. Gagie , G. Navarro, and Y. Nekrich . Ecient fullycompressed sequence representations . Algorithmica . To appear.

[2] G. Battaglia , D. Cangelosi , R. Grossi , and N. Pisanti . Masking patterns in sequences: A new class of motif discovery with don't cares . Theoretical

[7] M. David , M. Dzamba , D. Lister , L. Ilie , and M. Brudno . SHRiMP2: Sensitive yet practical short read mapping . Bioinformatics , 27 : 1011 - 1012 , 2011 .

[8] L. Egidi and G. Manzini . Better spaced seeds using quadratic residues . Journal of Compututer and System Sciences , 79 : 1144 - 1155 , 2013 .

[9] N. Homer , B. Merriman , and S. F. Nelson . BFAST: An alignment tool for large scale genome resequencing . PLOS One , 4 : e7767 , 2009 .

[10] P. Horton , S. M. Kielbasa , and M. C. Frith . DisLex: A tranformation for discontiguous sux array construction . In Proceedings of the Workshop on Knowledge, Language, and Learning in Bioinformatics (KLLBI) , pages 1 - 11 , 2008 .

[11] L. Ilie , S. Ilie , S. Khoshraftar , and A. Mansouri Bigvand . Seeds for e↵ective oligonucleotide design . BMC Genomics , 12 : 280 , 2011 .

[12] S. M. Kielbasa , R. Wan , K. Sato , P. Horton , and M. C. Frith . Adaptive seeds tame genomic sequence comparison . Genome Research , 21 : 487 - 493 , 2011 .

[13] G. Kucherov , L. No´e, and M. A. Roytberg . A unifying framework for seed sensitivity and its application to subset seeds . Journal of Bioinformatics and Computational Biology , 4 : 553 - 570 , 2006 .

[14] B. Langmeand and S. L. Salzberg . Fast gapped-read alignment with Bowtie 2 . Nature Methods , 9 : 357 - 359 , 2012 .

[15] B. Ma , J. Tromp, and M. Li . PatternHunter: faster and more sensitive homology search . Bioinformatics , 18 : 440 - 445 , 2002 .

[16] T. Onodera and T. Shibuya . An index structure for spaced seed search . In Proceedings of the 22nd International Symposium on Algorithms and Computation (ISAAC) , pages 764 - 772 , 2011 .

[17] P. Peterlongo , N. Pisanti , F. Boyer , and M.-F. Sagot . Lossless filter for finding long multiple approximate repetitions using a new data structure, the bifactor array . In Proceedings of the 12th Symposium on String Processing and Information Retrieval (SPIRE) , pages 179 - 190 , 2005 .

[18] L. M. S. Russo and G. Tischler . Succinct gapped sux arrays . In Proceedings of the 17th Symposium on String Processing and Information Retrieval (SPIRE) , pages 290 - 294 , 2011 .

[19] Y. Sun and J. Buhler . Designing multiple simultaneous seeds for DNA similarity search . Journal of Computational Biology , 12 : 847 - 861 , 2005 .

[20] K. J. Supowit . Decomposing a set of points into chains, with applications to permutation and circle graphs . Information Processing Letters , 21 : 249 - 252 , 1985 .