gsufsort: constructing suffix arrays, LCP arrays and BWTs for string collections

Background The construction of a suffix array for a collection of strings is a fundamental task in Bioinformatics and in many other applications that process strings. Related data structures, as the Longest Common Prefix array, the Burrows–Wheeler transform, and the document array, are often needed to accompany the suffix array to efficiently solve a wide variety of problems. While several algorithms have been proposed to construct the suffix array for a single string, less emphasis has been put on algorithms to construct suffix arrays for string collections. Result In this paper we introduce gsufsort, an open source software for constructing the suffix array and related data indexing structures for a string collection with N symbols in O(N) time. Our tool is written in ANSI/C and is based on the algorithm gSACA-K (Louza et al. in Theor Comput Sci 678:22–39, 2017), the fastest algorithm to construct suffix arrays for string collections. The tool supports large fasta, fastq and text files with multiple strings as input. Experiments have shown very good performance on different types of strings. Conclusions gsufsort is a fast, portable, and lightweight tool for constructing the suffix array and additional data structures for string collections.


Background
The suffix array ( SA ) [1] is one of the most important data structures in string processing. It enables efficient pattern searching in strings, as well as solving many other string problems [2][3][4]. More space-efficient solutions for such problems are possible by replacing the suffix array with an index based on the Burrows-Wheeler transform ( BWT ) [5]. Many applications require additional data structures-most commonly, the longest common prefix ( LCP ) [6] array and the document array ( DA ) [7]-on top of SA or BWT . These structures, possibly stored in compressed form, serve as a basis for building modern compact full-text indices, which allow to efficiently preprocess and query strings in compact space.
There are several internal memory algorithms designed for constructing the suffix array and additional data structures when the input consists of a single string [8,9]. While less emphasis has been put on specialized algorithms for string collections, in many applications the input is composed by many strings, and a common approach is concatenating all strings into a single one and using a standard construction algorithm. However, this approach may deteriorate either the theoretical bounds or the practical behavior of construction algorithms due to, respectively, the resulting alphabet size or unnecessary string comparisons [10][11][12].
Textual documents and webpages are examples of widespread large string collections. In Bioinformatics, important problems on collections of sequences may be solved rapidly with a small memory footprint using the aforementioned data structures, for example, finding suffix-prefix overlaps for sequence assembly [13] cDNA sequences [14], finding repeats [15] and sequence matching [16].
In this paper we present gsufsort, an open source tool that takes a string collection as input, constructs its (generalized) suffix array and additional data structures, like the BWT , the LCP array, and the DA , and writes them directly to disk. This way, applications that rely on such data structures may either read them from disk or may easily include gsufsort as a component. Large collections, with up to 2 64 − d − 2 total letters in d strings, may be handled provided that there is enough memory. This tool is an extension of previous results [10], with new implementations of procedures to obtain the BWT and the generalized suffix array ( GSA ) from SA during output to disk, and with the implementation of a lightweight alternative to compute DA. Implementation gsufsort is implemented in ANSI C and requires a single Make command to be compiled. It may receive a collection of strings in fasta, fastq or raw ASCII text formats and computes SA and related data structures, according to input parameters. gsufsort optionally supports gzipped input data using zlib 1 and kseq 2 libraries. Setting command-line arguments allows selecting which data structures are computed and written on disk, and which construction algorithm is used (see below). Additionally, a function for loading pre-constructed data structures from disk is also provided.
Given a collection of d strings T 1 , T 2 , . . . , T d from an alphabet � = [1, σ ] of ASCII symbols, having lengths n 1 , n 2 , . . . , n d , the strings are concatenated into a single string T [0, N − 1] = T 1 $T 2 $ · · · $T d $# using the same separator $ and an end-marker #, such that $ and # do not occur in any string T i , and # < $ < α for any other symbol α ∈ � . The total length of T is d i=1 (n i + 1) + 1 = N. Before giving details on gsufsort implementation, we briefly recall some data structures definitions. For a string S of length n let the suffix starting at position i be denoted S i , 0 ≤ i ≤ n − 1 . The suffix array SA of a string S of length n is an array with a permutation of [0, n − 1] that gives the suffixes of S in lexicographic order. The length of the longest common prefix of strings R and S is denoted by lcp(R, S) . The LCP array for S gives the lcp between consecutive suffixes in the order of SA , that is where each entry (a, b) represents the suffix T a b , with 1 ≤ a ≤ d and 0 ≤ b ≤ n a − 1.
gsufsort uses algorithm gSACA-K [10] to construct SA for the concatenated string T [0, N − 1] , which breaks ties between equal suffixes from different strings T i and T j by their ranks, namely i and j. gSACA-K can also compute LCP and DA during SA construction, such that LCP values do not exceed separator symbols. gSACA-K runs in O(N) time using O(σ ) working space.
The BWT is calculated during the output to disk according to its well-known relation to SA [3] The generalized suffix array ( GSA ) can be computed by gsufsort from SA and DA during the output to disk, using the identity

Results
We compared our tool and mkESA. mkESA [17] is a fast suffix array construction software designed for bioinformatics applications.
We ran both versions of our tool, gsufsort and gsufsort-light, to build arrays GSA and LCP , while mkESA 3 was run to build arrays SA and LCP for the concatenation of all strings (using the same symbol as separators). The experiments were conducted on a single core of a machine with GNU/Linux (Debian 8, kernel 3.16.0-4, 64 bits) with an Intel Xeon E5-2630 2.40-GHz, 384 GB RAM and 13 TB SATA storage. The sources were compiled by GNU GCC version 4.8.4 with option -O3.
The collections we used in our experiments are described in Table 1. They comprise real DNAs, real proteins, documents, random DNA and random protein, and differ by their alphabet size and also by the maximum (1) and average lcp , which offer an approximation for suffix sorting difficulty. The results are shown in Table 2. The data shows a clear time/memory tradeoff for DNA sequences, gsufsort being faster while using approximately 1.25 more memory, gsufsort-light using slightly less memory then mkESA but taking more time. On proteins, gsufsort-light is only marginally slower than gsufsort but faster than mkESA. The authors of mkESA reported a 32% gain on a large protein dataset using 16 threads [17], but larger lcp values seem not to favor mkESA when compared to gsufsort-light, which is 47.9% faster on proteins and 12.9% faster on DNA.
The memory ratio (bytes/N) of gsufsort and gsufsort-light is constant, 21 and 17 bytes per input symbol respectively, corresponding to the space of the input string T (N bytes) plus the space for arrays SA and LCP (8N bytes each) and, only for gsufsort, the space for DA (4N bytes).
We have also evaluated the performance of gsufsort, gsufsort-light and mkESA on collections of random DNA and random protein sequences. The collections have a growing number of 1MB sequences. The running time in seconds and the peak memory usage in GB are shown in Fig. 1 (logarithmic scale). Using random sequences reduces the variation due to lcp among Table 1 Collections gutenberg are ASCII books in English from Project Gutenberg (http://www.guten berg.org); random-dna was generated with even sampling probability on the standard 4 letter alphabet; random-protein was generated with even sampling probability on the IUPAC 25 letter alphabet  collections. We can see a perfectly steady behavior of mkESA. While still O(N), gsufsort displays a deviation due to larger constants.

Conclusions
We have introduced gsufsort, a fast, portable, and lightweight tool for constructing the suffix array and additional data structures for string collections. gsufsort may be used to pre-compute indexing structures and write them to disk, or may be included as a component in different applications. As an additional advantage, gsufsort is not restricted to biological sequences, as it can process collections of strings over ASCII alphabets.