External memory BWT and LCP computation for sequence collections with applications

Background Sequencing technologies produce larger and larger collections of biosequences that have to be stored in compressed indices supporting fast search operations. Many compressed indices are based on the Burrows–Wheeler Transform (BWT) and the longest common prefix (LCP) array. Because of the sheer size of the input it is important to build these data structures in external memory and time using in the best possible way the available RAM. Results We propose a space-efficient algorithm to compute the BWT and LCP array for a collection of sequences in the external or semi-external memory setting. Our algorithm splits the input collection into subcollections sufficiently small that it can compute their BWT in RAM using an optimal linear time algorithm. Next, it merges the partial BWTs in external or semi-external memory and in the process it also computes the LCP values. Our algorithm can be modified to output two additional arrays that, combined with the BWT and LCP array, provide simple, scan-based, external memory algorithms for three well known problems in bioinformatics: the computation of maximal repeats, the all pairs suffix–prefix overlaps, and the construction of succinct de Bruijn graphs. Conclusions We prove that our algorithm performs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(n\, \mathsf {maxlcp})$$\end{document}O(nmaxlcp) sequential I/Os, where n is the total length of the collection and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathsf {maxlcp}$$\end{document}maxlcp is the maximum LCP value. The experimental results show that our algorithm is only slightly slower than the state of the art for short sequences but it is up to 40 times faster for longer sequences or when the available RAM is at least equal to the size of the input.


Introduction
A fundamental problem in bioinformatics is the ability to efficiently search into the billions of DNA sequences produced by NGS studies. The Burrows Wheeler transform (BWT) [1] is a well known structure which is the starting point for the construction of compressed indices for collections of sequences [2]. The BWT is often complemented with the longest common prefix (LCP) array [3] since the latter makes it possible to efficiently emulate Suffix Tree algorithms [4,5]. The construction of such data structures is a challenging problem. Although the final outcome is a compressed index, construction algorithms can be memory hungry and the necessity of developing lightweight algorithms was recognized since the very beginning of the field [6][7][8]. In lightweight algorithms it is assumed that the input and the output fit in main memory but the amount of additional working memory is sublinear with the size of the input.
When even lightweight algorithms do not fit in RAM, one has to resort to external or semi-external memory construction algorithms (see [9][10][11][12][13][14] and references therein). In the external memory model [15] it is assumed that the main memory grows at most polylogarithmically with the size of the input. In the semi-external model the main memory can grow linearly with the size of the input but part of the working data has to reside on disk. In both models the complexity of the algorithms is usually measured in terms of disk I/Os, since data transfer is much slower than CPU operations.
The space efficient computation of the BWT in main memory for a single sequence is well studied, and remarkable advances have been recently obtained [16,17]. Unfortunately, for external memory computation the situation is less satisfactory. For collections of sequences, the first external memory algorithm is the BCR algorithm described in [18] that computes the multi-string BWT for a collection of total size n, performing a number of sequential I/Os proportional to nK, where K is the length of the longest sequence in the collection. This approach is clearly not competitive when the sequences have non-homogeneous lengths, and it is far from the theoretical optimal even for sequences of equal length. Nevertheless, the simplicity of the algorithm makes it very effective for collections of relatively short sequences, and it has become the reference tool for this problem. This approach was later extended [19] to compute also the LCP values with the same asymptotic number of I/ Os. When computing also the LCP values, or when the input strings have different lengths, the algorithm uses O(m) words of RAM, where m is the number of input sequences.
In this paper, we present a new space-efficient algorithm for the computation of the BWT and LCP array for a collection of sequences in external or semi-external memory. Our algorithm takes the amount of available RAM as an input parameter, and tries to make the best use of it by splitting the input into subcollections sufficiently small so that it can compute their BWT in internal memory using an optimal linear time algorithm. Next, it merges the partial BWTs in external or semiexternal memory and in the process it also computes the LCP values. Since the LCP values are computed in a non-standard order, the algorithm is completed by an external memory mergesort procedure that computes the final LCP array. We show that our algorithm performs a number of sequential I/Os between O(n avelcp) and O(n maxlcp) , where avelcp and maxlcp are respectively the average and the maximum longest common prefix of the input sequences. To our knowledge, the only other known external memory algorithm for computing the BWT and LCP arrays of a collection of sequences is bwt-lcp-em, recently proposed in [20] that performs O(n maxlcp) sequential I/Os and uses O(m + K ) words of RAM, where K is the fixed string length.
In "Related approaches" section we compare our approach with the ideas behind these previous works, and in "Experiments" section we compare their performance in practice. The experimental results show that BCR is the fastest algorithm for relatively short sequences while our algorithm is the fastest when the average LCP of the collection is relatively small compared to the lengths of the sequences. Both our algorithm and BCR appear to be faster than the available implementation of bwt-lcpem, which is however only a prototype implementation with some limitations on the admissible inputs.
Another contribution of the paper, which follows from our first result, is the design of simple external memory algorithms for three well known problems related to genomic sequences, namely: (i) the computation of maximal repeats [21,22], (ii) the computation of the all pairs suffix-prefix overlaps [23][24][25], and (iii) the construction of succinct de Bruijn graphs [26][27][28]. Our external memory algorithms for these problems are derived from known internal memory algorithms, but they process the input data in a single sequential scan. In addition, for the problem of computing the all pairs suffix-prefix, we go beyond the recent solutions that compute all the overlaps [24,25,29,30], and we compute only the overlaps above a certain length, still spending constant time per reported overlap. Our algorithms take as input the BWT and LCP array, together with two additional arrays that our BWT construction algorithm can compute without any asymptotic slowdown.
Since problems on genomic sequences often involve huge datasets, it is certainly important to provide efficient external memory algorithms for the three problems described above. To our knowledge, only for the all pairs suffix-prefix problem there exists an external memory algorithm, namely the algorithm [30, Algorithm 2] that computes all the overlaps given the BWT, LCP, and Generalized Suffix Array of the input collection.

Background
Let s[1, n] denote a string of length n over an alphabet of constant size σ . As usual, we assume s[n] is a special symbol (end-marker) not appearing elsewhere in s and lexicographically smaller than any other symbol. We write s[i, j] to denote the substring s is the empty string. Given two strings s 1 and s 2 we write s 1 s 2 ( s 1 ≺ s 2 ) to denote that s 1 is lexicographically (strictly) smaller than s 2 . We denote by LCP(s 1 , s 2 ) the length of the longest common prefix between s 1 and s 2 .
The suffix array sa[1, n] associated to s is the permutation of [1, n] giving the lexicographic order of s 's suffixes, that is, The longest common prefix array lcp[1, n + 1] is defined for i = 2, . . . , n by the lcp array stores the length of the longest common prefix (LCP) between lexicographically consecutive suffixes. For convenience we define lcp [ Let s 1 [1, n 1 ], . . . , s k [1, n k ] denote a collection of strings such that s 1 [n 1 ] = $ 1 , . . . , s k [n k ] = $ k , where $ 1 < . . . < $ k are k symbols not appearing elsewhere in s 1 , . . . , s k and smaller than any other symbol. Let sa 1···k [1, n] denote the suffix array of the concatenation s 1 · · · s k of total length n = k h=1 n h . The multi-string BWT [19,31] of s 1 , . . . , s k , denoted by bwt 1···k [1, n] , is defined as Essentially bwt 1···k is a permutation of the symbols in s 1 , . . . , s k such that the position in bwt 1···k of s i [j] is given by the lexicographic rank of its context s i [j + 1, n i ] (or Figure 1 shows an example with k = 2 .
Notice that for k = 1 , this is the usual Burrows-Wheeler Transform [1].
Given the suffix array sa 1···k [1, n] of the concatenation s 1 · · · s k , we consider the corresponding LCP array lcp 1···k [1, n] defined as in (1) (see again Fig. 1). Note that, for i = 2, . . . , n , lcp 1···k [i] gives the length of the longest common prefix between the contexts of bwt 1···k [i] and bwt 1···k [i − 1] . We stress that all practical implementations use a single $ symbol as end-marker for all strings to avoid alphabet explosion, but end-markers from Fig. 1 LCP array and BWT for s 1 = abcab$ 1 and s 2 = aabcabc$ 2 , and multi-string BWT and corresponding LCP array for the same strings. Column id shows, for each entry of bwt 12 = bc$ 2 cc$ 1 aaaabbb whether it comes from s 1 or s 2 different strings are then sorted as described, i.e., on the basis of the index of the strings they belong to.

Computing multi-string BWTs
The gSACA-K algorithm [32], based on algorithm SACA-K [33], computes the suffix array for a string collection. Given a collection of strings of total length n, gSACA-K computes the suffix array for their concatenation in O(n) time using (σ + 1) log n additional bits (in practice, only 2KB are used for ASCII alphabets). gSACA-K is time and space optimal for alphabets of constant size σ = O(1) . The multi-string bwt 1···k of s 1 , . . . , s k can be easily obtained from the suffix array as in (2). gSACA-K can also compute the lcp array lcp 1···k still in linear time using only the additional space for the lcp values.

Merging multi-string BWTs
The Gap algorithm [34], based on an earlier algorithm by Holt and McMillan [35], is a simple procedure for merging multi-string BWTs. In its original formulation the Gap algorithm can also merge LCP arrays, but in this paper we compute LCP values using a different approach more suitable for external memory execution. We describe here only the main idea behind Gap and refer the reader to [34] for further details. For simplicity in the following we assume we are merging k single-string BWTs bwt 1 = bwt(s 1 ), . . . , bwt k = bwt(s k ) ; the algorithm does not change in the general case where the inputs are multi-string BWTs. Computing bwt 1···k amounts to sorting the symbols of bwt 1 , . . . , bwt k according to the lexicographic order of their contexts, where the context of symbol bwt j [i] is s j [sa j [i], n j ] . By construction, the symbols in each bwt j are already sorted by context, hence to compute bwt 1···k we only need to merge bwt 1 , . . . , bwt k without changing the relative order of the symbols within the input sequences.
The Gap algorithm works in successive iterations. During the h-th iteration it computes a vector Z (h) specifying how the entries of bwt 1 , . . . , bwt k should be merged to have them sorted according to the first h symbols of their context. Formally, for j = 1, . . . , k the vector Z (h) contains |bwt j | copies of the value j arranged so that the following property holds. Property 1 is equivalent to state that we can logically partition Z (h) into b(h) + 1 blocks such that each block corresponds to the set of symbols in bwt 1···k , whose contexts are prefixed by the same length-h string. The context of any symbol in block Z (h) [ℓ j + 1, ℓ j+1 ] is lexicographically smaller than the context of the symbols in block Z (h) [ℓ k + 1, ℓ k+1 ] with k > j ; within each block, if j 1 < j 2 the symbols of bwt j 1 precede those of bwt j 2 . We keep explicit track of such blocks using a bit array B[1, n + 1] such that at the end of iteration h it is B[i] � = 0 if and only if a block of Z (h) starts at position i. By Property 1, when all entries in B are nonzero, Z (h) describes how the bwt j ( j = 1, . . . , k ) should be merged to get bwt 1···k .

Related approaches
The strategy used by Gap to build multi-string BWTs follows the idea, introduced by [35,36], of merging partial BWTs, i.e. BWTs of subsets of the input collection. Interestingly, both previous algorithms for computing the BWT and LCP in external memory [19,20] are also based on a merging strategy but instead of merging partial BWTs, they merge the arrays L 1 , L 2 , L 3 , …, where L i consists of the symbols which are at distance i from the end of their respective strings. The symbols inside each L i are sorted as usual by context. In the example of Fig. 1, we would have L 1 = bc (since s 1 ends with b$ 1 and s 2 ends with c$ 2 ), L 2 = ab , (since s 1 ends with ab$ 1 and s 2 ends with bc$ 2 ), L 3 = ca and so on. Note that in L 3 c precedes a since c 's context ab$ 1 is lexicographically smaller than a 's context bc$ 2 . Clearly, merging the arrays L i yields the desired multi-string BWT and the authors of [19,20] show how to compute also the LCP array. The algorithms in [19,20] differ in how the merging is done: [19] uses a refinement of a technique introduced in [9,10], while [20] uses a refinement of Holt and McMillan merging strategy [35,36].

The eGap algorithm
The eGap algorithm for computing the multi-string BWT and LCP array in external memory works in three phases. First it builds multi-string BWTs for sub-collections in internal memory, then it merges these BWTs in external memory and generates the LCP values. Finally, it sorts the LCP values in external memory.

Phase 1: BWT computation
Given a collection of sequences s 1 , s 2 , . . . , s k , we split it into sub-collections sufficiently small that we can compute the multi-string SA for each one of them in internal memory using the gSACA-K algorithm. After computing each SA, we compute the corresponding multi-string BWT and write it to disk in uncompressed form using one byte per character.

Phase 2: BWT merging and LCP computation
This part is based on the Gap algorithm previously described but suitably modified to work efficiently in external memory (or in semi-external memory if there are at least n bytes of RAM). In the following we assume that the input consists of k BWTs bwt 1 , . . . , bwt k of total length n over an alphabet of size σ . The input BWTs are read from disk and never moved to internal memory. The algorithm initially sets Z (0) = 1 n 1 2 n 2 . . . k n k and B = 10 n−1 1 . Since the context of every symbol is prefixed by the same length-0 string (the empty string), initially there is a single block containing all symbols. At iteration h the algorithm computes Z (h) from Z (h−1) as follows (see also Fig. 2). We define an array contains the number of occurrences of characters smaller than c in bwt 1···k . F partitions Z (h) into σ buckets, one for each symbol. Using Z (h−1) we scan the partially merged BWT, and whenever we encounter the BWT character c coming from bwt ℓ , with ℓ ∈ {1, . . . , k} , we store it in the next free position of bucket c in Z (h) ; note that c is not actually moved, instead we write ℓ in its corresponding position in Z (h) . In our implementation, instead of using distinct arrays Z (0) , Z (1) , . . . we only use two arrays Z old and Z new , that are kept on disk. At the beginning of iteration h it is Z old = Z (h−1) and Z new = Z (h−2) ; at the end Z new = Z (h) and the roles of the two files are swapped. While Z old is accessed sequentially, Z new is updated sequentially within each bucket, that is within each set of positions corresponding to a given character. Since the boundary of each bucket is known in advance we logically split the Z new file in buckets and write to each one sequentially.
eGap computes LCP values exploiting the bitvector B used by Gap to mark the beginning of blocks (see Eq. 3) within each Z (h) (for simplicity the computation of B is not shown in Fig. 2 has been set at iteration h − 1 : thus we output to a temporary file F h−2 the pair �i, h − 2� to record that lcp 1···k [i] = h − 2 , and we set B x [i] = 1 so no pair for position i will be produced in the following itera- the pairs are written in increasing order of their first component, since B and B x are scanned sequentially. These temporary files will be merged in Phase 3 to produce the LCP array. As proven in [34,Lemma 7], if at iteration h of the Gap algorithm we set B[i] = 1 , then at any iteration g ≥ h + 2 processing the entry Z (g) [i] will not change the arrays Z (g+1) and B. Since the roles of the Z old and Z new files are swapped at each iteration, and at iteration h we scan can be added to a set of irrelevant ranges that the algorithm may skip in successive iterations (irrelevant ranges are defined in terms of the array B x as opposed to the array B, since before skipping an irrelevant range we need to update both Z old and Z new ). We read from one file the ranges to be skipped at the current iteration and simultaneously write to another file the ranges to be skipped at the next iteration (note that irrelevant ranges are created and consumed sequentially). Since skipping a single irrelevant range takes O(k + σ ) time, an irrelevant range is stored only if its size is larger than a given threshold t and we merge consecutive irrelevant ranges whenever possible.
In our experiments we used t = max(256, k + σ ) . In the worst case the space for storing irrelevant ranges could be O(n) but in actual experiments it was always less than 0.1n bytes. As in the Gap algorithm, when all entries in B are nonzero, Z old describes how the BWTs bwt j ( j = 1, . . . , k ) should be merged to get bwt 1···k , and a final sequential scan of the input BWTs along with Z old allows to write bwt 1···k to disk, in sequential order. Our implementation can merge at most 2 7 = 128 BWTs at a time because we use 7 bits to store each entry of Z old and Z new . These arrays are maintained on disk in two separate files; the additional bit of each byte are used to keep the current and the next copy of B. The bit array B x is stored separately in a file of size n/8 bytes. To merge a set of k > 128 we split the input in subsets of cardinality 128 and merge them in successive rounds. In practice, the algorithm merges the multi-string BWTs produced by Phase 1. In our experiments the maximum number of sub-collections was 21.
Semi-external version We have also implemented a semi-external version of the merge algorithm that uses n bytes of RAM. The i-th byte is used to store . This version can sort at most 2 3 = 8 BWTs simultaneously; to sort k BWTs it performs log 8 k merging rounds. Although performing more rounds is clearly more expensive, this version stores in RAM all the arrays that are modified and reads from disk only the input BWTs and is therefore significantly faster.

Analysis
During Phase 1, gSACA-K computes the suffix array for a sub-collection of total length m using 9m bytes (8m bytes for sa and 1m bytes for the text). If the available RAM is M, the input is split into subcollections of size ≈ M/9 . Since gSACA-K runs in linear time, if the input collection has total size n, Phase 1 takes O(n) time overall.
A single iteration of Phase 2 consists of a complete scan of Z (h−1) except for the irrelevant ranges. Since the algorithm requires maxlcp iterations, without skipping the irrelevant ranges the algorithm would require maxlcp sequential scans of O(n) items. Reasoning as in [34,Theorem 8] we get that by skipping irrelevant ranges the overall amount of data directly read/written by the algorithm is O(n avelcp) items where avelcp is the arithmetic average of the entries in the final LCP array. However, if we reason in terms of disk blocks, every time we skip an irrelevant range we discard the current block and load a new one (unless the beginning of the new relevant range is inside the same block; in that case or if the beginning of the new relevant range is in the block immediately following the current one, skipping the irrelevant range does not save any I/O). We can upper bound this extra cost, with an overhead of O(1) blocks for each irrelevant range skipped. Summing up, if the total number of skipped ranges is Ir and each disk block consists of B words, the I/O complexity of Phase 2 according to the theoretical model in [15] is O(Ir + n avelcp/(B log n)) block I/Os (under the reasonable assumptions that the alphabet is constant, each entry in Z takes constant space, and we need a constant number of merge rounds). Although the experiments in "Experiments" section suggest that in practice Ir is small, for simplicity and uniformity with the previous literature we upper bound the cost of Phase 2 with O(n maxlcp) sequential I/Os (corresponding to O(n maxlcp/(B log n)) block I/Os).
Phase 3 takes O(⌈log maxlcp⌉) rounds; each round merges LCP files by sequentially reading and writing O(n) bytes of data. The overall cost of Phase 3 is therefore O(n log maxlcp) sequential I/Os. In our experiments we used = 256 ; since in our tests maxlcp < 2 16 two merging rounds were always sufficient.

Experiments
In this section we report on an experimental study comparing between the eGap algorithm and the other known external memory tools computing the BWT and LCP arrays of sequence collections. We implemented eGap in ANSI C based on the code of Gap [34] and gSACA-K [32]. eGap source code is freely available at https ://githu b.com/felip elouz a/egap/. All tested algorithms were compiled with GNU GCC ver. Memory setting To make a realistic external memory experimental setting one has to use an amount of RAM smaller than the size of the data. Indeed, if more RAM is available, even if the algorithm is supposedly not using it, the operating system will use it to temporarily store disk data and the algorithm will be no longer really working in external memory. This phenomenon will be apparent also from our experiments. For these reasons we reduced the available RAM to simulate three different scenarios: (i) input data 4 times larger than the available RAM, (ii) input data of approximately the same size as the RAM, and (iii) input data 4 times smaller than the RAM. We evaluated these scenarios with the complete 8 GB datasets from Table 1 (with  2 GB, 8 GB, and 32 GB RAM), and with the datasets trimmed to 1 GB (hence with 256 MB, 1 GB, and 4 GB RAM). The RAM was limited at boot time to a value equal to the amount assigned to the algorithm plus a small extra amount for the operating system (14 MB for the 256 MB instance and 64 MB for the others).

Comparison with the existing algorithms
We compared eGap with the algorithm BCR [19] which is the current state of the art for BWT/LCP computation for collections of sequences. We used the bcr-lcp implementation from [38] since the previous implementation mentioned in [19] did not compute the LCP values correctly. We tested also the recently proposed algorithm bwt-lcp-em [20] using the code from [39]. As a reference we also tested the algorithm eGSA [14] using the code from [40]. eGSA computes the Suffix and LCP Arrays for collections of sequences in external memory: the disadvantage of this algorithm is that working with the Suffix Array could involve transferring to/from disk a much larger amount of data.
Limitations We tested bwt-lcp-em only on the short 1 GB dataset since the implementation in [39] only supports collections of at most 2 GB and with strings of at most 253 symbols. We tested eGSA only with memory scenario (iii) (input data 4 times smaller than the RAM) since it was already observed in [14] that eGSA 's running time degrades when the RAM is restricted to the input size. Finally, we could not test bcr-lcp on the pacbio 1 GB dataset since it stopped with an internal error after four days of computation. This is probably due to the presence of very long strings in the dataset since bcr-lcp was originally conceived for collections of short/medium length strings. The corresponding entries are marked as "failed" in Fig. 3. For the larger 8 GB datasets we stopped the experiments that did not complete after six days of CPU time, corresponding to 60 microseconds per input symbol. The corresponding entries are marked with " > 60 " in Fig. 3. Note that both bwt-lcpem and bcr-lcp are active projects, so some of the limitations reported here could have been solved after our experiments were completed.
Results: The results of our experiments are summarized in Fig. 3. The bar plots on the left are for the 1 GB datasets showing the running time as function of the available RAM; the diagrams on the right are for the 8 GB datasets. The results show that for memory scenarios (i) and (ii) eGap and bcr-lcp have the better performance, whereas for scenario (iii) eGap and eGSA are the best options. The performance of bwt-lcp-em improves with the RAM size, but it is still 12 times slower than eGap for the short datasets with 4 GB of RAM.
The above results are in good accordance with the theoretical analysis. bcr-lcp complexity is O(n maxlen) sequential I/Os while eGap and bwt-lcp-em both take     O(n maxlcp) sequential I/Os. For the short and long datasets the maximum length and the maximum LCP coincide and we see that when the available memory is only one fourth of the input size bcr-lcp is clearly the fastest option: indeed it is up to a factor 2.6 faster than eGap. This is no longer true when the available memory is equal or larger than the input size: in this case eGap is the fastest, probably because of its ability to exploit all the available memory using a semi-external strategy whenever possible. When the available memory is larger than the input size or for the pacbio.1000 dataset which has a very large maxlen then eGap is up to 40 times faster than bcr-lcp. Note that, in accordance with our heuristic analysis, eGap 's running time per input byte appears to be roughly proportional to the average LCP of the collection. If we look at the datasets pacbio and pacbio.1000 we see that they have widely different maximum LCPs, yet their running times are very close similarly to their average LCPs. Note that in the scenario (iii) eGSA is often the fastest algorithm and its running time appears to be less influenced by the size of the average or maximum LCP. Another advantage is that it also computes the Suffix Array, but it has the drawback of using a large amount of disk working space: 340 GB for a 8 GB input vs 56 GB used by eGap.
We conclude that, although eGap is not always the fastest algorithm, its running time is never too far from that of the best algorithm. In addition, eGap is the only algorithm that was able to complete all computations in all memory models. Although it was devised as an external memory algorithm, its ability to switch to a semiexternal strategy if the memory is available makes it a very flexible tool. The comparison with the other algorithms in this setting is indeed not completely fair, since none of them is designed to take the available memory as a parameter in order to make the best use of it. Note that, as the available memory increases, all algorithms become faster because the operating system uses the RAM as a buffer but the speed improvement is different for different algorithms.

Relative performance of eGap's building blocks
We evaluated the percentage of time spent by each phase of eGap and their efficiency (percentage the CPU was busy) on the 8 GB datasets in the memory scenarios considered above, thus with RAM limited to (i) 1 GB, (ii) 8 GB, and (iii) 32 GB.
The results in Fig. 4 show that Phase 2 of eGap dominates the algorithm in general. The second phase took about 95% , 85% and 50% of the total time in scenarios (i), (ii), and (iii) respectively. If we look at the efficiency of the single phases, we see that they all improve with the RAM size. However, we notice that for any given memory scenario the efficiency of Phases 1 and 3 was almost the same for the different datasets, while Phase 2 has a different behavior. For the short and long datasets with 8 GB and 32 GB RAM, we see that Phase 2 efficiency is very close to Phase 1's, while there is a sharp drop when using 2 GB RAM. For the pacbio datasets, the drop in Phase 2 efficiency is significant already when we use 8 GB RAM.

Applications
In this section we show that the eGap algorithm, in addition to the BWT and LCP arrays, can output additional information useful to design efficient external memory algorithms for three well known problems on sequence collections: (i) the computation of maximal repeats, (ii) the all pairs suffix-prefix overlaps, and (iii) the construction of succinct de Bruijn graphs. For these problems we describe algorithms which are derived from known (internal memory) algorithms suitably modified so that they process the input data in a single sequential scan.
Our first observation is that eGap can also output the array which provides, for each bwt entry, the id of the sequence to which that entry belongs. In information retrieval this is usually called the Document Array, so in the following we will denote it by da . In Phase 1 the gSACA-K algorithm can compute the da together with the lcp and bwt using only additional 4n bytes of space to store the da entries. These partial da 's can be merged in Phase 2 using the Z new array in the same way as the BWT entries. In the following we use bwt , lcp , and da to denote the multistring BWT, LCP and Document Array of a collection of m sequences of total length n. We write s to denote the concatenation s 1 · · · s m and sa to denote the suffix array of s . We will use s and sa to describe and prove the correctness of our algorithms, but neither s nor sa are used in the computations.

Computation of maximal repeats
Different notions of maximal repeats have been used in the bioinformatics literature to model different notions of repetitive structure (see for example [21,22]). We use a notion of maximal repeat from [41,Chap. 7]: we say that a string α is a Type 1 maximal repeat if α occurs in the collection at least twice and every extension, i.e. cα or αc with c ∈ , occurs fewer times. We consider also a more restrictive notion: we say that a string α is a Type 2 maximal repeat if α occurs in the collection at least twice and every extension of α occurs at most once.
To compute Type 1 maximal repeats the crucial observation is that there is a substring of length ℓ that prefixes sa entries j, j + 1, . . . , i (and no others) iff lcp[k] ≥ ℓ for k = j + 1, . . . , i , and both lcp[j] and lcp[i + 1] are smaller . Thus, to ensure that we have found a Type 1 maximal repeat we only need to check that bwt[j − 1, i] contains at least two distinct characters. To efficiently check this latter condition, for each stack entry �j, lcp[h]� we maintain a bit vector b j of size σ keeping track of the distinct characters in the array bwt from position j − 1 to the next stack entry, or to the last seen position for the entry at the top of the stack. When �j, lcp[h]� is popped from the stack its bit vector is or-ed to the previous stack entry in constant time; if �j, lcp[h]� is popped from the stack and immediately replaced with �j, lcp[i]� its bit vector survives as it is (essentially because it is associated with an index, not with a stack entry). Clearly, maintaining the bit vector does not increase the asymptotic cost of the algorithm.
Since at each step we insert at most one entry on the stack, the overall cost of our algorithm is O(n) time. The algorithm uses a stack of size bounded by O(maxlcp) words. For most applications maxlcp ≪ n so it should be feasible to keep the stack in RAM. However, since a stack can also be implemented in external memory in O(1) amortized time per operation [42], we can state the following result. , n] are prefixed by the same string α of length lcp[j + 1] and every extension αc occurs at most once. If this is the case, then α is a Type 2 maximal repeat if all characters in bwt[j, i] are distinct since this ensures that also every extension cα occurs at most once. In order to detect this situation, as we scan the lcp array we maintain a candidate pair �j + 1, lcp[j + 1]� such that j + 1 is the largest index seen so far for which lcp[j] < lcp[j + 1] . When we establish a candidate at j + 1 as above, we initialize to zero a bit vector b of size σ setting to 1 only entries bwt[j] and bwt[j + 1] . As long as the following values lcp[j + 2], lcp[j + 3], . . . are equal to lcp[j + 1] we go on updating b and if the same position is marked twice we discard �j + 1, lcp[j + 1]� . If we reach an index i + 1 such that lcp[i + 1] > lcp[j + 1] , we update the candidate to �i + 1, lcp[i + 1]� and reinitialize b. If we reach i + 1 such that lcp[i + 1] < lcp[j + 1] and �j + 1, lcp[j + 1]� has not been discarded, then a repeat of Type 2 (with i − j + 1 repetitions) has been located. Note that when our algorithms discover Type 1 or Type 2 maximal repeats we know the repeat length and the number of occurrences so one can easily filter out noninteresting repeats (too short or too frequent). In some applications, for example the MUMmer tool [43], one is interested in repeats that occur in at least r distinct input sequences, maybe exactly once for each sequence. Since for these applications the number of input sequences is relatively small, we can handle these requirements by simply scanning the da array simultaneously with the lcp and bwt arrays and keeping track of the sequences associated to a maximal repeat using a bit vector (or a unionfind structure) as we do with characters in the bwt.

All pairs suffix-prefix overlaps
In this problem we want to compute, for each pair of sequences s i s j , the longest overlap between a suffix of s i and a prefix of s j . Our solution is inspired by the algorithm in [24] which in turn was derived by an earlier Suffix-tree based algorithm [23]. The algorithm in [24] solves the problem using a Generalized Enhanced Suffix array (consisting of the arrays sa , lcp , and da ) in O(n + m 2 ) time, which is optimal since n is the size of the input and there are m 2 longest overlaps. However, for large collections it is natural to consider the problem of reporting only the overlaps larger than a given threshold τ still spending O(n) time plus constant time per reported overlap. Our algorithm solves this more challenging problem.
In the following we say that a suffix starting at sa[i] is special iff it is a prefix of the suffix starting at sa[i + 1] , not considering the end-marker. This is equivalent to state that s[sa[i] + lcp[i + 1]] = $ . For example, in Fig. 1 (right) the special suffixes are ab$ 1 , abc$ 2 , abcab$ 1 b$ 1 , bc$ 2 , bcab$ 1 , c$ 2 , cab$ 1 . Notice that a special suffix starting at sa[i] has the form v$ with |v| = lcp[i + 1] ; clearly only if sa[i] is special then v can be a suffix-prefix overlap. Note also that any suffix $ is always trivially special.
To efficiently solve the suffix-prefix overlaps problem, we modify Phase 2 of our algorithm so that it outputs also the bit array xlcp such that xlcp[i] = 1 iff the suffix starting at sa[i] is special. To this end, we maintain an additional length-n bit array S such that, at the end of iteration h, S[i] = 1 if and only if the suffix starting at sa[i] is special and it has length less than h, again not considering the end-marker symbol. The array S is initialized at the end of iteration h = 1 as S = 1 k 0 n−k , consistently with the fact that in the final suffix array the first k contexts are strings consisting of just an end-marker, that are special suffixes and the only suffixes of length 0.
During iteration h, we update S as follows. With reference to the code in Fig. 2, whenever we use entry Step 5 of Fig. 2, then the context corresponding to position j is cv$ . Since the context contains the end-marker, j is the correct lexicographic position of cv$ which is therefore the suffix corresponding to sa The computation of the overlaps is done as in [24]. When the scan reaches position i, we check whether bwt[i] = $ . If this is the case, then s[sa[i], n] is prefixed by the whole sequence s da [i] , hence the longest overlap between a prefix of s da[i] and a suffix of s k is given by the element currently at the top of stack[k] , since by construction these stacks only contain special suffixes whose overlap with s[sa[i], n] is larger than τ . Note that using lcpStack and top we can directly access the stacks whose top element corresponds to an overlap with s da[i] larger than τ , hence the time spent in this phase is proportional to the number of reported overlaps. As in [24] some care is required to handle the case in which the whole string s da[i] is a suffix of another sequence, but this can be done without increasing the overall complexity as in [24]. Since we spend constant time for reported overlap and amortized constant time for scanned position the overall cost of the algorithm, in addition to the scanning of the bwt /lcp/xlcp/da arrays, is O(n + E τ ) , where E τ is the number of suffix-prefix overlaps greater than τ . Since all stacks can be implemented in external memory spending amortized constant time per operation, we only need to store in RAM top and daPtr that overall take O(m + maxlcp) words.

Theorem 3
Our algorithm computes all suffix-prefix overlaps longer than τ in time O(n + E τ ) , where E τ is the number of reported overlaps, using O(m + maxlcp) words of RAM and executing a single scan of the arrays bwt , lcp , da and xlcp .

Construction of succinct de Bruijn graphs
A recent remarkable application of compressed data structures is the design of efficiently navigable succinct representations of de Bruijn graphs [26][27][28]. Formally, a de Bruijn graph for a collection of strings consists of a set of vertices representing the distinct k-mers appearing in the collection, with a directed edge (u, v) iff there exists a (k + 1)-mer α in the collection such that α[1, k] is the k-mer associated to u and α[2, k + 1] is the k-mer associated to v.
The starting point of all de Bruijn graphs succinct representation is the BOSS representation [28], so called from the authors' initials. For simplicity we now describe the BOSS representation of a k-order de Bruijn graph using the lexicographic order of k-mers, instead of the co-lexicographic order as in [28], which means we are building the graph with the direction of the arcs reversed. This is not a limitation since arcs can be traversed in both directions (or we can apply our construction to the input sequences reversed).
Consider the N k-mers appearing in the collection sorted in lexicographic order. For each k-mer α i consider the array C i of distinct characters c ∈ ∪ {$} such that cα i appears in the collection. The concatenation W = C 1 C 2 · · · C N is the first component of the BOSS representation. The second component is a binary array last , with |last| = |W | , such that last[j] = 1 iff W [j] is the last entry of some array C i . Clearly, there is a bijection between entries in W and graph edges; in the array last each sequence 0 i 1 ( i ≥ 0 ) corresponds to the outgoing edges of a single vertex with outdegree i + 1 . Finally, the third component is a binary array W − , with |W − | = |W | , such that W − [j] = 1 iff W [j] comes from the array C i , where α i is the lexicographically smallest k-mer prefixed by α i [1, k − 1] and preceded by W[j] in the collection. This means that α i is the lexicographically smallest k-mer with an outgoing edge reaching the node associated to kmer W [j]α i [1, k − 1] . Note that the number of 1 's in last and W − is exactly N, i.e. the number of nodes in the de Bruijn graph.
We now show how to compute W , last and W − by a sequential scan of the bwt and lcp array. The crucial observation is that the suffix array range prefixed by the same k-mer α i is identified by a range If, in addition to the bwt and lcp arrays, we also scan the da array, then we can keep track of which sequences contain any given graph edge and therefore obtain a succinct representation of the colored de Bruijn graph [44]. Finally, we observe that if our only objective is to build the k-order de Bruijn graph, then we can stop the phase 2 of our algorithm after the k-th iteration. Indeed, we do not need to compute the exact values of LCP entries greater than k, and also we do not need the exact BWT but only the BWT characters sorted by their length k context.

Conclusions
In this paper we have described eGap, a new algorithm for the computation of the BWT and LCP arrays of large collection of sequences. Depending on the amount of available memory, eGap uses an external or semi-external strategy for computing the BWT and LCP values. An experimental comparison of the available tools for BWT and LCP arrays computation shows that eGap is the fastest tool in many scenarios and was the only tool capable of completing the computation within a reasonable time frame for all kind of input data.
Another important feature of eGap is that, in addition to the BWT and LCP array, it can compute, without any asymptotic slowdown, two additional arrays that provide important information about the substrings of the input collection. We show how to use such information to design efficient external memory algorithms for three important problems for biosequences, namely the computation of maximal repeats, the computation of the all pairs suffix-prefix overlaps, and the construction of succinct de Bruijn graphs. Overall our results confirm the importance of the BWT and LCP arrays beyond their use for the construction of compressed full text indexes. This is in accordance with other recent results that have shown of they can be used directly to discover structural information on the underlying collection (see [45][46][47] and references therein).