Relative Suffix Trees

Abstract Suffix trees are one of the most versatile data structures in stringology, with many applications in bioinformatics. Their main drawback is their size, which can be tens of times larger than the input sequence. Much effort has been put into reducing the space usage, leading ultimately to compressed suffix trees. These compressed data structures can efficiently simulate the suffix tree, while using space proportional to a compressed representation of the sequence. In this work, we take a new approach to compressed suffix trees for repetitive sequence collections, such as collections of individual genomes. We compress the suffix trees of individual sequences relative to the suffix tree of a reference sequence. These relative data structures provide competitive time/space trade-offs, being almost as small as the smallest compressed suffix trees for repetitive collections, and competitive in time with the largest and fastest compressed suffix trees.


INTRODUCTION
The suffix tree [1] is one of the most powerful bioinformatic tools to answer complex queries on DNA and protein sequences [2,3,4].A serious problem that hampers its wider use on large genome sequences is its size, which may be 10-20 bytes per character.In addition, the non-local access patterns required by most interesting problems solved with suffix trees complicate secondary-memory deployments.This problem has led to numerous efforts to reduce the size of suffix trees by representing them using compressed data structures [5,6,7,8,9,10,11,12,13,14,15,16,17], leading to compressed suffix trees (CST).Currently, the smallest CST is the so-called fully-compressed suffix tree (FCST) [10,14], which uses 5 bits per character (bpc) for DNA sequences, but takes milliseconds to simulate suffix tree navigation operations.In the other extreme, Sadakane's CST [5,11] uses about 12 bpc and operates in microseconds, and even nanoseconds for the simplest operations.
A space usage of 12 bpc may seem reasonable to handle, for example, one human genome, which has about 3.1 billion bases: it can be operated within a RAM of 4.5 GB (the representation contains the sequence as well).However, as the price of sequencing has fallen, sequencing the genomes of a large number of individuals has become a routine activity.The 1000 Genomes Project [18] sequenced the genomes of several thousand humans, while newer projects can be orders of magnitude larger.This has made the development of techniques for storing and analyzing huge amounts of sequence data flourish.
Just storing 1000 human genomes using a 12 bpc CST requires almost 4.5 TB, which is much more than the amount of memory available in a commodity server.Assuming that a single server has 256 GB of memory, we would need a cluster of 18 servers to handle such a collection of CSTs (compared to over 100 with classical suffix tree implementations!).With the smaller (and much slower) FCST, this would drop to 7-8 servers.It is clear that further space reductions in the representation of compressed suffix trees would lead to reductions in hardware, communication, and energy costs when implementing complex searches over large genomic databases.
An important characteristic of those large genome databases is that they usually consist of the genomes of individuals of the same or closely related species.This implies that the collections are highly repetitive, that is, each genome can be obtained by concatenating a relatively small number of substrings of other genomes and adding a few new characters.When repetitiveness is considered, much higher compression rates can be obtained in compressed suffix trees.For example, it is possible to reduce the space to 1-2 bpc (albeit with operation times in the milliseconds) [13], or to 2-3 bpc with operation times in the microseconds [15].Using 2 bpc, our 1000 genomes could be handled with just 3 servers with 256 GB of memory.We note, however, that these CSTs index the whole collection and not individual sequences, which makes a difference in the types of queries that can be answered.This also makes a distributed implementation less obviously scalable.
Compression algorithms best capture repetitiveness by using grammar-based compression or Lempel-Ziv compression. 7In the first case [19,20] one finds a context-free grammar that generates (only) the text collection.The more repetitive the collection is, the smaller the grammar becomes.Rather than compressing the text directly, the current CSTs for repetitive collections [13,15] apply grammar-based compression on the data structures that simulate the suffix tree.Grammar-based compression yields relatively easy direct access to the compressed sequence [21], which makes it attractive compared to Lempel-Ziv compression [22], despite the latter generally using less space.
Lempel-Ziv compression cuts the collection into phrases, each of which has already appeared earlier in the collection.To extract the content of a phrase, one may have to recursively extract the content at that earlier position, following a possibly long chain of indirections.So far, the indexes built on Lempel-Ziv compression [23] or on combinations of Lempel-Ziv and grammar-based compression [24,25,26] support only pattern matching, which is just one of the wide range of functionalities offered by suffix trees.The high cost to access the data at random positions lies at the heart of the research on indexes built on Lempel-Ziv compression.
A simple way out of this limitation is the socalled relative Lempel-Ziv (RLZ) compression [27], where one of the sequences is represented in plain form and the others can only take phrases from that reference sequence.This enables immediate access for the symbols inside any copied phrase (as no transitive referencing exists) and, at least if a good reference sequence has been found, offers compression competitive with the classical Lempel-Ziv.In our case, taking any random genome per species as the reference is good enough; more sophisticated techniques have been studied [28,29,30].Structures for direct access [31,32] and even for pattern matching [33,34] have been developed on top of RLZ.
In this paper, we develop a CST by augmenting the relative FM-index [34] with structures based on RLZ.On a collection of human genomes, we achieve less than 3 bpc and operate within microseconds.This performance is comparable to that of a previous CST for this scenario [15], but our CSTs have a different functionality.
We have a separate CST for each sequence, instead of a single CST for their concatenation.Depending on the application, one kind of CST or the other is necessary.
Our compressed suffix tree, called RST, follows a trend of CSTs [6,7,9,8,11,13] that use only a patternmatching index (called suffix array) and an array with the length of the longest common prefix between each suffix and the previous one in lexicographic order (called LCP).We use the relative FM-index as our suffix array, and compress LCP using RLZ.On top of the RLZ phrases we build a tree of range minima that enables fast range minimum queries, as well as next-and previous-smallervalue queries, on LCP [13].All the CST functionality is built on those queries [6].Our main algorithmic contribution is this RLZ-based representation of the LCP array with the required extra functionality.
Another approach to compressing a repetitive collection while supporting interesting queries is to build an automaton that accepts the sequences in the collection, and then index the state diagram as an directed acyclic graph (DAG); see, for example, [35,36,37] for recent discussions.The first data structure to take this approach was the Generalized Compressed Suffix Array (GCSA) [38,37], which was designed for pangenomics so queries can return information about sequences not in the collection but that can be obtained from those in the collection by recombination.The FMindex of an alignment (FMA) [39,40] is similar to the GCSA but indexes only the sequences in the collection: whereas the GCSA conceptually embeds the automaton in a de Bruijn graph, the FMA embeds it in a coloured de Bruijn graph [41], preserving its specificity.Both the GCSA and FMA are practical but neither support the full functionality of a suffix tree.The precursor to the FMA, the suffix tree of an alignment (STA) [42,43], allows certain disjunctions in the suffix tree's edge labels in order to reduce the size of the tree while maintaining its functionality.It differs from our approach because it is a single tree for the whole collection and not a separate one for each sequence; also, unlike the FMA, the STA has not been implemented.Both the STA and FMA divide the sequences in the collection into regions of variation and conserved regions, and depend on the conserved regions being long enough that they can be distinguished from each other and the variations.This dependency makes these structures vulnerable to even a small change in even one sequence to an otherwiseconserved region, which could hamper their scalability.
For indexing purposes, we often consider text strings T [1, n] that are terminated by an endmarker T [n] = $ = 0 not occurring elsewhere in the text.Binary sequences are sequences over the alphabet {0, For any binary sequence B[1, n], we define the subsequence S[B] of string S [1, n] as the concatenation of the characters s i with B Contiguous subsequences S[i, j] are called substrings.Substrings of the form S [1, j] and S[i, n], i, j ∈ [1, n], are called prefixes and suffixes, respectively.We define the lexicographic order among strings in the usual way.

Full-text indexes
The suffix tree (ST) [1] of text T is a trie containing the suffixes of T , with unary paths compacted into single edges.Because the degree of every internal node is at least two, there can be at most 2n − 1 nodes, and the suffix tree can be stored in O(n log n) bits.In practice, this is at least 10n bytes for small texts [44], and more for large texts as the pointers grow larger.If v is a node of a suffix tree, we write π(v) to denote the concatenation of the labels of the path from the root to v.
Suffix arrays (SA) [45] were introduced as a spaceefficient alternative to suffix trees.The suffix array SA T [1, n] of text T is an array of pointers to the suffixes of the text in lexicographic order. 8In its basic form, the suffix array requires n log n bits in addition to the text, but its functionality is more limited than that of the suffix tree.In addition to the suffix array, many algorithms also use the inverse suffix array ISA [1, n], with SA[ISA[i]] = i for all i.
Let lcp(S 1 , S 2 ) be the length of the longest common prefix (LCP) of strings S 1 and S 2 .
The LCP array [45] LCP [1, n] of text T stores the LCP lengths for lexicographically adjacent suffixes of T as Abouelhoda et al. [46] showed how traversals on the suffix tree could be simulated using the suffix array, the LCP array, and a representation of the suffix tree topology based on lcp-intervals, paving the way for more space-efficient suffix tree representations.

Compressed text indexes
Data structures supporting rank and select queries over sequences are the main building blocks of compressed text indexes.If S is a sequence, we define rank c (S, i) as the number of occurrences of character c in the prefix S [1, i], while select c (S, j) is the position of the occurrence of rank j in sequence S. A bitvector is a representation of a binary sequence supporting fast rank and select queries.Wavelet trees (WT) [47] use bitvectors to support rank and select on general sequences.
The Burrows-Wheeler transform (BWT) [48] is a reversible permutation BWT[1, n] of text T .It is defined as ). Originally intended for data compression, the Burrows-Wheeler transform has been widely used in space-efficient text indexes, because it shares the combinatorial structure of the suffix tree and the suffix array.
Let LF be a function such that where C[c] is the number of occurrences of characters with lexicographical values smaller than c in BWT.The inverse function of LF is Ψ, with ), where c is the largest character value with C[c] < i.With functions Ψ and LF, we can move forward and backward in the text, while maintaining the lexicographic rank of the current suffix.If the sequence S is not evident from the context, we write LF S and Ψ S .
Compressed suffix arrays (CSA) [49,50,51] are text indexes supporting a functionality similar to the suffix array.This includes the following queries: i) find(P ) = [sp, ep] determines the lexicographic range of suffixes starting with pattern P [1, ]; ii) locate(sp, ep) = SA[sp, ep] returns the starting positions of these suffixes; and iii) extract(i, j) = T [i, j] extracts substrings of the text.In practice, the find performance of CSAs can be competitive with suffix arrays, while locate queries are orders of magnitude slower [52].Typical index sizes are less than the size of the uncompressed text.
The FM-index (FMI) [50] is a common type of compressed suffix array.A typical implementation [53] stores the BWT in a wavelet tree [47].The index implements find queries via backward searching.Let [sp, ep] be the lexicographic range of the suffixes of the text starting with suffix P [i + 1, ] of the pattern.We can find the range matching suffix P [i, ] with a generalization of function LF as We support locate queries by sampling some suffix array pointers.If we want to determine a value SA[i] that has not been sampled, we can compute it as SA[i] = SA[j] + k, where SA[j] is a sampled pointer found by iterating LF k times, starting from position i.Given sample interval d, the samples can be chosen in suffix order, sampling SA[i] at positions divisible by d, or in text order, sampling T [i] at positions divisible by d and marking the sampled SA positions in a bitvector.Suffix-order sampling requires less space, often resulting in better time/space trade-offs in practice, while text-order sampling guarantees better worst-case performance.We also sample the ISA pointers for extract queries.To extract T [i, j], we find the nearest sampled pointer after T [j], and traverse backwards to T [i] with function LF.
Compressed suffix trees (CST) [5] are compressed text indexes supporting the full functionality of a suffix tree (see Table 1).They combine a compressed suffix array, a compressed representation of the LCP array, and a compressed representation of suffix tree topology.For the LCP array, there are several common representations: • LCP-byte [46] stores the LCP array as a byte array.
If LCP[i] < 255, the LCP value is stored in the byte array.Larger values are marked with a 255 in the byte array and stored separately.As many texts produce small LCP values, LCP-byte usually requires n to 1.5n bytes of space.• We can store the LCP array by using variablelength codes.LCP-dac uses directly addressable codes [54] for the purpose, resulting in a structure that is typically somewhat smaller and somewhat slower than LCP-byte.
the LCP array stored in text order and used as the array can be stored as a bitvector of length 2n in 2n+o(n) bits.If the text is repetitive, run-length encoding can be used to compress the bitvector to take even less space [6].Because accessing PLCP uses locate, it is much slower than the above two encodings.
Suffix tree topology representations are the main difference between the various CST proposals.While the compressed suffix arrays and the LCP arrays are interchangeable, the tree representation determines how various suffix tree operations are implemented.There are three main families of compressed suffix trees: • Sadakane's compressed suffix tree (CST-Sada) [5] uses a balanced parentheses representation for the tree.Each node is encoded as an opening parenthesis, followed by the encodings of its children and a closing parenthesis.This can be encoded as a bitvector of length 2n , where n is the number of nodes, requiring up to 4n + o(n) bits.CST-Sada tends to be larger and faster than the other compressed suffix trees [11,13].• The fully compressed suffix tree (FCST) of Russo et al. [10,14] [7,9,8,11,13], the CST-NPR is perhaps the most practical compressed suffix tree.
For typical texts and component choices, the size of compressed suffix trees ranges from the 1.5n to 3n bytes of CST-Sada to the 0.5n to n bytes of FCST [11,13].There are also some CST variants for repetitive texts, such as versioned document collections and collections of individual genomes.Abeliuk et al. [13] developed a variant of CST-NPR that can sometimes be smaller than n bits, while achieving performance similar to the FCST.Navarro and Ordóñez [15] used grammar-based compression for the tree representation of CST-Sada.The resulting compressed suffix tree (GCT) requires slightly more space than the CST-NPR of Abeliuk et al., while being closer to the non-repetitive CST-Sada and CST-NPR in performance.

Relative Lempel-Ziv
Relative Lempel-Ziv (RLZ) parsing [27] compresses target sequence S relative to reference sequence R.The target sequence is represented as a concatenation of z phrases w i = (p i , i , c i ), where p i is the starting position of the phrase in the reference, i is the length of the copied substring, and c i is the mismatch character.If phrase w i starts from position p in the target, then S[p , p The shortest RLZ parsing of the target sequence can be found in (essentially) linear time.The algorithm builds a CSA for the reverse of the reference sequence, and then parses the target sequence greedily by using backward searching.If the edit distance between the reference and the target is s, we need at most s phrases to represent the target sequence.On the other hand, because the relative order of the phrases can be different in sequences R and S, the edit distance can be much larger than the number of phrases in the shortest RLZ parsing.
In a straightforward implementation, the phrase pointers p i and the mismatch characters c i can be stored in arrays W p and W c .These arrays take z log|R| and z log σ bits, respectively.To support random access to the target sequence, we can encode phrase lengths as a bitvector W of length |S| [27]: we set W is the first character of a phrase.The bitvector requires z log n z + O(z) bits if we use the sdarray representation [55].To extract S[j], we first determine the phrase w i , with i = rank 1 (W , j).If W [j + 1] = 1, we return the mismatch character W c [i]. Otherwise we determine the phrase offset with a select query, and return the character Ferrada et al. [32] showed how, by using relative TABLE 1.Typical compressed suffix tree operations.

Root()
The root of the tree.

Leaf(v)
Is node v a leaf?Ancestor(v, w) Is node v an ancestor of node w?

Count(v)
Number of leaves in the subtree with v as the root.

Locate(v)
Pointer to the suffix corresponding to leaf v.

Parent(v)
The parent of node v.

FChild(v)
The first child of node v in alphabetic order.

NSibling(v)
The next sibling of node v in alphabetic order.LCA(v, w) The lowest common ancestor of nodes v and w.

SDepth(v)
String depth: Length = |π(v)| of the label from the root to node v.

TDepth(v)
Tree depth: The depth of node v in the suffix tree.
The highest ancestor of node v with string depth at least d.
The ancestor of node v with tree depth d.
Suffix link iterated k times.

Child(v, c)
The child of node v with edge label starting with character c.
pointers instead of absolute pointers, we can avoid the use of select queries.They also achieved better compression of DNA collections, in which most of the differences between the target sequences and the reference sequence are single-character substitutions.
If most of the differences are single-character substitutions, p i+1 will often be p i + i +1.This corresponds to W r [i+1] = W r [i] with relative pointers, making run-length encoding of the pointer array worthwhile.
When we sort the suffixes in lexicographic order, substitutions in the text move suffixes around, creating insertions and deletions in the suffix array and related structures.In the LCP array, an insertion or deletion affecting LCP[i] can also change the value of LCP[i + 1].Hence RLZ with relative pointers is not enough to compress the LCP array.

Cox et al. [56] modified Ferrada et al.'s version of RLZ
to handle other small variations in addition to singlecharacter substitutions.After adding a phrase to the parse, we look ahead a bounded number of positions to find potential phrases with a relative pointer W r [i] close to the previous explicit relative pointer W r [j].If we can find a sufficiently long phrase this way, we encode the pointer differentially as W r [i] − W r [j].Otherwise we store W r [i] explicitly.We can then save space by storing the differential pointers separately using less bits per pointer.Because there can be multiple mismatch characters between phrases i and i + 1, we also need a prefix-sum data structure L for finding the range W c [a, b] containing the mismatches.Cox et al. showed that their approach compresses both DNA sequences and LCP arrays better than Ferrada et al.'s version, albeit with slightly slower random access.We refer the reader to their paper for more details of their implementation.

RELATIVE FM-INDEX
The relative FM-index (RFM) [34] is a compressed suffix array of a sequence relative to the CSA of another sequence.The index is based on approximating the longest common subsequence (LCS) of BWT R and BWT S , where R is the reference sequence and S is the target sequence, and storing several structures based on the common subsequence.Given a representation of BWT R supporting rank and select, we can use the relative index RFM S|R to simulate rank and select on BWT S .

Basic index
Assume that we have found a long common subsequence of sequences X and Y .
We call positions X[i] and Y [j] lcs-positions, if they are in the common subsequence.If B X and B Y are the binary sequences marking the common subsequence ( ), we can move between lcs-positions in the two sequences with rank and select operations.
In its most basic form, the relative FM-index RFM S|R only supports find queries by simulating rank queries on BWT S .It does this by storing Align(BWT R , BWT S ) and the complements LCS(BWT R ) and LCS(BWT S ) of the common subsequence.The lcs-bitvectors are compressed using entropy-based compression [57], while the complements are stored in structures similar to the reference BWT R .
To compute rank c (BWT S , i), we first determine the number of lcs-positions in BWT S up to position S[i] with k = rank 1 (B BWT S , i).Then we find the lcs-position k in BWT R with j = select 1 (B BWT R , k).With these positions, we can compute

Relative select
We can implement the entire functionality of a compressed suffix array with rank queries on the BWT.However, if we use the CSA in a compressed suffix tree, we also need select queries to support forward searching with Ψ and Child queries.We can always implement select queries by binary searching with rank queries, but the result will be much slower than the rank queries.
A faster alternative to support select queries in the relative FM-index is to build a relative select structure rselect [58].Let F X be a sequence consisting of the characters of sequence X in sorted order.Alternatively, as well as the C array C LCS for the common subsequence.
To compute select c (BWT S , i), we first determine how many of the first i occurrences of character c are lcs-positions with whether the occurrence we are looking for is an lcsposition or not.If it is, we find the position in ), and then map j to select c (BWT S , i) by using Align(BWT R , BWT S ).Otherwise we find the occurrence in LCS(BWT S ) with j = select c (LCS(BWT S ), i − k), and return select c (BWT S , i) = select 0 (B BWT S , j).

Full functionality
If we want the relative FM-index to support locate and extract queries, we cannot build it from any common subsequence of BWT R and BWT S .We need a bwtinvariant subsequence [34], where the alignment of the BWTs is also an alignment of the original sequences.Definition 3.1.Let X be a common subsequence of BWT R and BWT S , and let BWT R [i R ] and BWT S [i S ] be the lcs-positions corresponding to X for all positions i, j ∈ {1, . . ., |X|}.
In addition to the structures already mentioned, the full relative FM-index has another pair of lcs-bitvectors, Align(R, S), which marks the bwt-invariant subsequence in the original sequences.
To compute the answer to a locate(i) query, we start by iterating BWT S backwards with LF queries, until we find an lcs-position BWT S [i ] after k steps.Then we map position i to the corresponding position j in BWT R by using Align(BWT R , BWT S ).Finally we determine SA R [j ] with a locate query in the reference index, and map the result to SA S [i ] by using Align(R, S). 10 The result of the locate(i) query is The ISA S [i] access required for extract queries is supported in a similar way.We find the lcs-position S[i + k] for the smallest k ≥ 0, and map it to the corresponding position R[j] by using Align(R, S).Then we determine ISA R [j + 1] by using the reference index, and map it back to ISA S [i + k + 1] with Align(BWT R , BWT S ).Finally we iterate BWT S k + 1 steps backward with LF queries to find ISA S [i].
If the target sequence contains long insertions not present in the reference, we may also want to include some SA and ISA samples for querying those regions.

Finding a bwt-invariant subsequence
With the basic relative FM-index, we approximate the longest common subsequence of BWT R and BWT S by partitioning the BWTs according to lexicographic contexts, finding the longest common subsequence for each pair of substrings in the partitioning, and concatenating the results.The algorithm is fast, easy to parallelize, and quite space-efficient.As such, RFM construction is practical, having been tested with datasets of hundreds of gigabytes in size.
To find a bwt-invariant subsequence, we start by matching each suffix of the reference sequence with the lexicographically nearest suffixes of the target sequence.In contrast to the original algorithm [34], we only match suffixes that are lexicographically adjacent in the mutual suffix array of the two sequences.
Instead of using the mutual suffix array, we can use CSA R , CSA S , and the merging bitvector B R,S of length |RS|.We set B R,S [i] = 1, if SA RS [i] points to a suffix of S. We can build the merging bitvector in O(|S| • t LF ) time, where t LF is the time required for an LF query, by extracting S from CSA S and backward searching for it in CSA R [59] Our next step is building the match arrays left and right, which correspond to the arrays A[•] [2] and A[•] [1] in the original algorithm.This is done by traversing CSA R backwards from ISA R [|R|] = 1 with LF queries and following the left and the right matches of the current suffix.During the traversal, we maintain the invariant j = SA R [i] with (i, j) ← (LF R (i), j − 1).If suffix R[j, |R|] has a left (right) match, we use the shorthand l(j) = rank 1 (B R,S , select 0 (B R,S , i) − 1) (r(j) = rank 1 (B R,S , select 0 (B R,S , i) + 1)) to refer to its position in CSA S .
We say that suffixes R[j, |R|] and R[j + 1, |R|] have the same left match if l(j) = LF S (l(j + 1)).Let R[j, |R|] to R[j + , |R|] be a maximal run of suffixes having the same left match, with suffixes R[j, |R|] to R[j + − 1, |R|] starting with the same characters as their left matches. 11We find the left match of suffix R[j, |R|] as j = SA S [l(j)] by using CSA S , and set left[j, j + − 1] = [j , j + − 1].The right match array right is built in a similar way.
The match arrays require 2|R| log|S| bits of space.If sequences R and S are similar, the runs in the arrays tend to be long.Hence we can run-length encode the match arrays to save space.The traversal takes O(|R| • (t LF + t rank + t select ) + rd • t LF ) time, where t rank and t select denote the time required by rank and select operations, r is the number of runs in the two arrays, and d is the suffix array sample interval in CSA S . 12 The final step is finding the longest increasing subsequence X of arrays left and right, which corresponds to a common subsequence of R and S.More precisely, we want to find a binary sequence B R [1, |R|], which marks the common subsequence in R, and an integer sequence X, which contains the positions of the common subsequence in S. The goal is to make sequence X strictly increasing and as long as possible, with X[rank 1 (B R , i)] being either left[i] or right[i].This can be done in O(|R| log|R|) time with O(|R| log|R|) bits of additional working space with a straightforward modification of the dynamic programming algorithm for finding the longest increasing subsequence.While the dynamic programming tables can be run-length encoded, the time and space savings are negligible or even non-existent in practice.
As sequence X is strictly increasing, we can convert it into binary sequence B S [1, |S|], marking the values in sequence X with 1-bits.Afterwards, we can consider the binary sequences B R and B S as the lcs-bitvectors Align(R, S).Because every suffix of R starts with the same character as its matches stored in the left and right 11 The first character of a suffix can be determined by using the C array. 12 The time bound assumes text-order sampling.
arrays, subsequences R[B R ] and S[B S ] are identical.As each suffix R[i, |R|] with B R [i] = 1 is paired with its left match or right match in sequence S, no other suffix of R or S is lexicographically between the two paired suffixes.For any i, let i R = select 1 (B R , i) and i S = select 1 (B S , i) be the lcs-positions of rank i.Then, for 1 ≤ i, j ≤ |X|, which is equivalent to the condition in Definition 3.1.We can convert Align(R, S) to Align(BWT R , BWT S ) in O((|R| + |S|) • t LF ) time by traversing CSA R and CSA S backwards.The resulting subsequence of BWT R and BWT S is bwt-invariant.
Note that the full relative FM-index is more limited than the basic index, because it does not handle substring moves very well.Let R = xy and S = yx, for two random sequences x and y of length n/2 each.Because BWT R and BWT S are very similar, we can expect to find a common subsequence of length almost n.On the other hand, the length of the longest bwtinvariant subsequence is around n/2, because we can either match the suffixes of x or the suffixes of y in R and S, but not both.

RELATIVE SUFFIX TREE
The relative suffix tree (RST) is a CST-NPR of the target sequence relative to a CST of the reference sequence.It consists of two major components: the relative FM-index with full functionality and the relative LCP (RLCP) array.The optional relative select structure can be generated or loaded from disk to speed up algorithms based on forward searching.The RLCP array is based on RLZ parsing, while the support for nsv/psv/rmq queries is based on a minima tree over the phrases.
We make a similar observation in the relative setting.If target sequence S is similar to the reference sequence R, then their LCP arrays should also be similar.If there are long identical ranges LCP R [i, i + k] = LCP S [j, j + k], the corresponding DLCP ranges DLCP R [i + 1, i + k] and DLCP S [j + 1, j + k] are also identical.Hence we can use RLZ parsing to compress either the original LCP array or the DLCP array.
While the identical ranges are a bit longer in the LCP array, we opt to compress the DLCP array, because it behaves better when there are long repetitions in the sequences.In particular, assembled genomes often have long runs of character N , which correspond to regions of very large LCP values.If the runs are longer in the target sequence than in the reference sequence, the RLZ parsing of the LCP array will have many mismatch characters.The corresponding ranges in the DLCP array typically consist of values {−1, 0, 1}, making them much easier to compress.
We create an RLZ parsing of DLCP S relative to DLCP R , while using LCP R as the reference afterwards.The reference is stored in a structure we call slarray, which is a variant of LCP-byte.[46].Small values LCP R [i] < 255 are stored in a byte array, while large values LCP R [i] ≥ 255 are marked with a 255 in the byte array and stored separately.To quickly find the large values, we also build a rank 255 structure over the byte array.The slarray provides reasonably fast random access and very fast sequential access to the underlying array.
The RLZ parsing produces a sequence of phrases w i = (p i , i , c i ) (see Section 2.3; since we are using Cox et al.'s version, c i is now a string).Because some queries involve decompressing an entire phrase, we limit the maximum phrase length to 1024.Phrase lengths are encoded in the W bitvector in the usual way.We convert the mismatching DLCP values c i into absolute LCP values in the mismatch array W c , and store it as an slarray.The mismatch values are used as absolute samples for the differential encoding.
To access LCP S [j], we determine the phrase w i as usual, and check whether we should return a mismatch character.If so, we compute which one using a prefix sum query on L, and return it.If not, we determine the starting positions r i and s i of the phrase w i in the reference and the target, respectively.We can then compute the solution as To speed up the computation, we enforce that each RLZ phrase ends with at least one mismatch character, so LCP S [s i − 1] is readily available.After finding LCP S [j], accessing LCP S [j − 1] and LCP S [j + 1] is fast, as long as we do not cross phrase boundaries.

Supporting nsv/psv/rmq queries
Suffix tree topology can be inferred from the LCP array with range minimum queries (rmq) and next/previous smaller value (nsv/psv) queries [6].
Some suffix tree operations are more efficient if we also support next/previous smaller or equal value (nsev/psev) queries [13].Query nsev(i) (psev(i)) finds the next (previous) value smaller than or equal to LCP [i].
In order to support the queries, we build a 64ary minima tree over the phrases of the RLZ parsing.Each leaf node stores the smallest LCP value in the corresponding phrase, while each internal node stores the smallest value in the subtree.Internal nodes are created and stored in a levelwise fashion, so that each internal node, except perhaps the rightmost one of each level, has 64 children.
We encode the minima tree as two arrays.The smallest LCP values are stored in M LCP , which we encode as an slarray.Plain array M L stores the starting offset of each level in M LCP , with the leaves stored starting from offset M L [1] = 1.If i is a minima tree node located at level j, the corresponding minimum value is M LCP [i], the parent of the node is A range minimum query rmq(sp, ep) starts by finding the minimal range of phrases w l , . . ., w r covering the query and the maximal range of phrases w l , . . ., w r contained in the query (note that l ≤ l ≤ l + 1 and r − 1 ≤ r ≤ r).We then use the minima tree to find the leftmost minimum value and find the leftmost occurrence LCP[i] = j in phrase w k .If l < l and M LCP [l] ≤ j, we decompress phrase w l and find the leftmost minimum value LCP[i ] = j (with i ≥ sp) in the phrase.If j ≤ j, we update (i, j) ← (i , j ).Finally we check phrase w r in a similar way, if r > r and M LCP [r] < j.The answer to the range minimum query is LCP[i] = j, so we return (i, j). 13inally, the particular case where no phrase is contained in [sp, ep] is handled by sequentially scanning one or two phrases in LCP.
The remaining queries are all similar to each other.In order to answer query nsv(i), we start by finding the phrase w k containing position i, and then determining LCP[i].Next we scan the rest of the phrase to see whether there is a smaller value LCP[j] < LCP[i] later in the phrase.If so, we return (j, LCP[j]).Otherwise we traverse the minima tree to find the smallest k > k with M LCP [k ] < LCP[i].We decompress phrase w k , find the leftmost position j with LCP[j] < LCP[i], and return (j, LCP[j]).

EXPERIMENTS
We have implemented the relative suffix tree in C++, extending the old relative FM-index implementation. 14he implementation is based on the Succinct Data Structure Library (SDSL) 2.0 [60].Some parts of the implementation have been parallelized using OpenMP and the libstdc++ parallel mode.
As our reference CSA, we used the succinct suffix array (SSA) [53,61] implemented using SDSL components.Our implementation is very similar to csa wt in SDSL, but we needed better access to the internals than what the SDSL interface provides.SSA encodes the Burrows-Wheeler transform as a Huffmanshaped wavelet tree, combining fast queries with size close to the order-0 empirical entropy.This makes it the index of choice for DNA sequences [52].In addition to the plain SSA with uncompressed bitvectors, we also used SSA-RRR with entropy-compressed bitvectors [57] to highlight the the time-space trade-offs achieved with better compression We sampled SA in suffix order and ISA in text order.In SSA, the sample intervals were 17 for SA and 64 for ISA.In RFM, we used sample interval 257 for SA and 512 for ISA to handle the regions that do not exist in the reference.The sample intervals for suffix order sampling were primes due to the long runs of character N in the assembled genomes.If the number of long runs of character N in the indexed sequence is even, the lexicographic ranks of almost all suffixes in half of the runs are odd, and those runs are almost completely unsampled.This can be avoided by making the sample interval and the number of runs relatively prime.
The experiments were done on a system with two 16-core AMD Opteron 6378 processors and 256 GB of memory.The system was running Ubuntu 12.04 with Linux kernel 3.2.0.We compiled all code with g++ version 4.9.2.We allowed index construction to use multiple threads, while confining the query benchmarks to a single thread.As AMD Opteron uses a nonuniform memory access architecture, accessing memory controlled by the same physical CPU is faster than accessing remote memory controlled by another CPU.In order to ensure that all data structures are in local memory, we set the CPU affinity of the query benchmarks with the taskset utility.
As our target sequence, we used the maternal haplotypes of the 1000 Genomes Project individual NA12878 [62].As the reference sequence, we used the 1000 Genomes Project version of the GRCh37 assembly of the human reference genome. 15Because NA12878 is female, we also created a reference sequence without chromosome Y.
In the following, a basic FM-index is an index supporting only find queries, while a full index also supports locate and extract queries.

Indexes and their sizes
Table 2 lists the resource requirements for building the relative indexes, assuming that we have already built the corresponding non-relative structures for the sequences.As a comparison, building an FM-index for a human genome typically takes 16-17 minutes and 25-26 GB of memory.While the construction of the basic RFM index is highly optimized, the other construction algorithms are just the first implementations.Building the optional rselect structures takes 4 minutes using two threads and around 730 megabytes (|R| + |S| bits) of working space in addition to RFM and rselect.
The sizes of the final indexes are listed in Table 3.The full RFM is over twice the size of the basic index, but still 3.3--3.7 times smaller than the full SSA-RRR and 4.6--5.3times smaller than the full SSA.The RLCP array is 2.7 times larger than the RFM index with the full human reference and 1.5 times larger with the female reference.Hence having a separate female reference is worthwhile, if there are more than a few female genomes among the target sequences.The optional rselect structure is almost as large as the basic RFM index.
Table 4 lists the sizes of the individual components of the relative FM-index.Including the chromosome Y in the reference increases the sizes of almost all relative components, with the exception of LCS(BWT S ) and Align(R, S).In the first case, the common subsequence still covers approximately the same positions in BWT S as before.In the second case, chromosome Y appears in bitvector B R as a long run of 0-bits, which compresses well.The components of a full RFM index are larger than the corresponding components of a basic RFM index, because the bwt-invariant subsequence is shorter than the approximate longest common subsequence (see Table 2).
The size breakdown of the RLCP array can be seen in Table 5. Phrase pointers and phrase lengths take space proportional to the number of phrases.As there are more mismatches between the copied substrings with the full human reference than with the female reference, the absolute LCP values take a larger proportion of the total space with the full reference.Shorter phrase length increases the likelihood that the minimal LCP value in a phrase is a large value, increasing the size of the minima tree.

Query times
Average query times for the basic operations can be seen in Tables 6 and 7.The results for LF and Ψ queries in the full FM-indexes are similar to the earlier ones with basic indexes [58].Random access to the RLCP array is about 30 times slower than to the LCP array, while sequential access is 10 times slower.The nsv, psv, and rmq queries are comparable to 1-2 random accesses to the RLCP array.
We also tested the locate performance of the full RFM index, and compared it to SSA and SSA-RRR.We built the indexes with SA sample intervals 7, 17, 31, 61, and 127, using the reference without chromosome Y for RFM. 16The ISA sample interval was the maximum of 64 and the SA sample interval.We extracted 2 million random patterns of length 32, consisting of characters ACGT , from the target sequence, and measured the total time taken by find and locate queries.The results can be seen in Figure 1.While SSA and SSA-RRR 16 With RFM, the sample intervals apply to the reference SSA.

Synthetic collections
In order to determine how the differences between the reference sequence and the target sequence affect the size of relative structures, we built RST for various synthetic datasets.We took a 20 MB prefix of the human reference genome as the reference sequence, and generated 25 target sequences with every mutation rate p ∈ {0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1}.A total of 90% of the mutations were single-character substitutions, while 5% were insertions and another 5% deletions.The length of an insertion or deletion was k ≥ 1 with probability 0.2 • 0.8 k−1 .The results can be seen in Figure 2 (left).The size of the RLCP array grew quickly with increasing mutation rates, peaking at p = 0.01.At that point, the average length of an RLZ phrase was comparable to what could be found in the DLCP arrays of unrelated DNA sequences.With even higher mutation rates, the phrases became slightly longer due to the smaller average LCP values.The RFM index, on the other hand, remained small until p = 0.003.Afterwards, the index started growing quickly, eventually overtaking the RLCP array.
We also compared the size of the relative suffix tree to GCT [15], which is essentially a CST-Sada for repetitive collections.While the structures are intended for different purposes, the comparison shows how much additional space is used for providing access to the suffix trees of individual datasets.We chose to skip the CST-NPR for repetitive collections [13], as its implementation was not stable enough.
Figure 2 (right) shows the sizes of the compressed suffix trees.The numbers for RST include individual indexes for each of the 25 target sequences as well as the reference data, while the numbers for GCT are for a single index containing the 25 sequences.With low mutation rates, RST was not much larger than GCT.The size of RST starts growing quickly at around p = 0.001, while the size of GCT stabilizes at 3-4 bpc.

Suffix tree operations
In the final set of experiments, we compared the performance of RST to the SDSL implementations of various compressed suffix trees.We used the maternal haplotypes of NA12878 as the target sequence and the human reference genome without chromosome Y as the reference sequence.We built RST, CST-Sada, CST-NPR, and FCST for the target sequence.CST-Sada uses Sadakane's compressed suffix array (CSA-Sada) [49] as its CSA, while the other SDSL implementations use SSA.We used PLCP as the LCP encoding with both CST-Sada and CST-NPR, and also built CST-NPR with LCP-dac.
We used three algorithms for the performance comparison.
The first algorithm is preorder traversal of the suffix tree using SDSL iterators (cst dfs const forward iterator).
The iterators use operations Root, Leaf, Parent, FChild, and NSibling, though Parent queries are rare, as the iterators cache the most recent parent nodes.
The other two algorithms find the maximal substrings of the query string occurring in the indexed text, and report the lexicographic range for each such substring.This is a key task in common problems such as We used the paternal haplotypes of chromosome 1 of NA12878 as the query string in the maximal substrings algorithms.Because some tree operations in the SDSL compressed suffix trees take time proportional to the depth of the current node, we truncated the runs of character N in the query string into a single character.Otherwise searching in the deep subtrees would have made some SDSL suffix trees much slower than RST.
The results can be seen in Table 8.RST was 1.8 times smaller than FCST and several times smaller than the other compressed suffix trees.In depth-first traversal, RST was 4 times slower than CST-NPR and about 15 times slower than CST-Sada.FCST was orders of magnitude slower, managing to traverse only 5.3% of the tree before the run was terminated after 24 hours.
It should be noted that the memory access patterns of traversing CST-Sada, CST-NPR, and RST are highly local.Traversal times are mostly based on the amount of computation done, while memory latency is less important than in the individual query benchmarks.In RST, the algorithm is essentially the following: 1) compute rmq in the current range; 2) proceed recursively to the left subinterval; and 3) proceed to the right subinterval.This involves plenty of redundant work, as can be seen by comparing the traversal time (0.90 μs per node) to sequential RLCP access (0.017 μs per position).A faster algorithm would decompress large parts of the LCP array at once, build the corresponding subtrees in postorder [46], and traverse the resulting trees.
RST with rselect is as fast as CST-Sada in the forward algorithm, 1.8-2.7 times slower than CST-NPR, and 4.1 times faster than FCST.Without the additional structure, RST becomes 2.6 times slower.As expected [64], the backward algorithm is much faster than the forward algorithm.CST-Sada and RST, which combine slow backward searching with a fast tree, have similar performance to FCST, which combines fast searching with a slow tree.CST-NPR is about an order of magnitude faster than the others in the backward algorithm.

DISCUSSION
We have introduced relative suffix trees (RST), a new kind of compressed suffix tree for repetitive sequence collections.Our RST compresses the suffix tree of an individual sequence relative to the suffix tree of a reference sequence.It combines an already known relative suffix array with a novel relative-compressed longest common prefix representation (RLCP).When the sequences are similar enough (e.g., two human genomes), the RST requires about 3 bits per symbol on each target sequence.This is close to the space used by the most space-efficient compressed suffix trees designed to store repetitive collections in a single tree, but the RST provides a different functionality as it indexes each sequence individually.The RST supports query and navigation operations within a few microseconds, which is competitive with the largest and fastest compressed suffix trees.
While our RST implementation provides competitive time/space trade-offs, there is still much room for Most importantly, some of the construction algorithms require significant amounts of time and memory.In many places, we have chosen simple and fast implementation options, even though there could be alternatives that require significantly less space without being too much slower.
Our RST is a relative version of the CST-NPR.Another alternative for future work is a relative CST-Sada, using RLZ compressed bitvectors for suffix tree topology and PLCP.q q q q q q SSA SSA−RRR RFM q q q q q q q q RFM RLCP Reference RST an internal node of the suffix tree, = |π(v)| the string depth of node v, and SA[sp, ep] the corresponding suffix array interval.The following properties hold for the lcp-interval LCP[sp, ep]

Definition 3 . 2 .
Let R and S be two sequences, and let SA = SA RS and ISA = ISA RS .The left match of

FIGURE 1 .
FIGURE 1.Average find and locate times in microseconds per occurrence for 2 million patterns of length 32 with a total of 255 million occurrences on NA12878 relative to the human reference genome without chromosome Y.Left: Query time vs. suffix array sample interval.Right: Query time vs. index size in bits per character.

TABLE 2 .
Sequence lengths and resources used by index construction for NA12878 relative to the human reference genome with and without chromosome Y.Approx and Inv denote the approximate LCS and the bwt-invariant subsequence.Sequence lengths are in millions of base pairs, while construction resources are in minutes of wall clock time and gigabytes of memory.

TABLE 3 .
Various indexes for NA12878 relative to the human reference genome with and without chromosome Y.The total for RST includes the full RFM.Index sizes are in megabytes and in bits per character.RFM used 5.4-7.6 microseconds per occurrence more than SSA, resulting in slower growth in query times.In particular, RFM with reference sample interval 31 was faster than SSA with sample interval 61.

TABLE 4 .
Breakdown of component sizes in the RFM index for NA12878 relative to the human reference genome with and without chromosome Y in bits per character.

TABLE 5 .
Breakdown of component sizes in the RLCP array for NA12878 relative to the human reference genome with and without chromosome Y.The number of phrases, average phrase length, and the component sizes in bits per character."Parse" contains Wr and W , "Literals" contains Wc and L, and "Tree" contains M LCP and ML.

TABLE 7 .
Query times in microseconds in the LCP array (slarray) and the RLCP array for NA12878 relative to the human reference genome with and without chromosome Y.For the random queries, the query times are averages over 100 million queries.The range lengths for the rmq queries were 16 k (for k ≥ 1) with probability 0.5 k .For sequential access, we list the average time per position for scanning the entire array.

TABLE 8 .
Index size in bits per character vs. mutation rate for 25 synthetic sequences relative to a 20 MB reference.Compressed suffix trees for the maternal haplotypes of NA12878 relative to the human reference genome without chromosome Y.Component choices; index size in bits per character; average time in microseconds per node for preorder traversal; and average time in microseconds per character for finding maximal substrings shared with the paternal haplotypes of chromosome 1 of NA12878 using forward and backward algorithms.The figures in parentheses are estimates based on the progress made in the first 24 hours.