Indexing labeled sequences

Background Labels are a way to add some information on a text, such as functional annotations such as genes on a DNA sequences. V(D)J recombinations are DNA recombinations involving two or three short genes in lymphocytes. Sequencing this short region (500 bp or less) produces labeled sequences and brings insight in the lymphocyte repertoire for onco-hematology or immunology studies. Methods We present two indexes for a text with non-overlapping labels. They store the text in a Burrows–Wheeler transform (BWT) and a compressed label sequence in a Wavelet Tree. The label sequence is taken in the order of the text (TL-index) or in the order of the BWT (TLBW-index). Both indexes need a space related to the entropy of the labeled text. Results These indexes allow efficient text–label queries to count and find labeled patterns. The TLBW-index has an overhead on simple label queries but is very efficient on combined pattern–label queries. We implemented the indexes in C++ and compared them against a baseline solution on pseudo-random as well as on V(D)J labeled texts. Discussion New indexes such as the ones we proposed improve the way we index and query labeled texts as, for instance, lymphocyte repertoire for hematological and immunological studies.


INTRODUCTION
Labels are a way to add some information on a text, as the semantics of words on an English sentence or functional annotations such as genes on a DNA sequences. Can we build an index that saves a labeled text like ACGCC : : : TTGA (of size 96), which have the label L 1 in the positions 12-41 and the label L 2 in the positions 56-96? We consider here the case where a same label can be given to different (but similar) patterns and can occur several times in the text. We introduce two indexes which store a labeled text and answers to position-label association queries. Those indexes share some ideas with the RL-FMI (Mäkinen & Navarro, 2004) which uses a Burrows-Wheeler transform (BWT) and a Wavelet Tree (WT). Using a somewhat similar organization, we index a labeled text. The following sections present the TL-and TL BW -indexes (text-label indexes) and their associated queries. The last section presents experimental results on simulated and real genomic data.
Let T = t 0 t 1 : : : t n -1 be a text of length n over an alphabet of size s. The text may be composed of several sequences, each sequence ending with the symbol $. Let L = {L 0 , L 1 : : : L l -1 } be a set of labels. A labeled text (T, A) is here a text with non-overlapping labels: a letter should have at most one label. Each position i of the text is labeled by exactly one label a i ∈ L ∪ {3}, where the special 3 label is put on every letter without label. A = a 0 a 1 : : : a n -1 is the label string. Figure 1 shows the text which will be used all over this article.
The BWT (Burrows & Wheeler, 1994) is a reversible algorithm which reorganizes the letters of a text. The transformed text, BWT(T), is the concatenation of the last letters of the lexicographically sorted rotations of the text. BWT(T) is easier to compress and can be stored using nH k (T) + o(n) bits, where H k is the kth order empirical entropy of the text T. The FM-index (Ferragina & Manzini, 2000) uses the BWT, a table C, where C [a] gives the number of letters lexicographically smaller than a, and a function Occ(a, i) giving the number of occurrences of a in BWT(T)[0, i]. The FM-index allows to efficiently search a pattern using backward search.
A WT (Grossi, Gupta & Vitter, 2003) is a binary tree storing a text, where each symbol from the alphabet corresponds to a leaf. The root is a bit vector where every position corresponds to the element it has to index. Any position marked 0 (respectively 1) corresponds to an element whose leaf is on the left (respectively right) descendant of the node. The process is repeated recursively until the leaves. For a text A of length a in an alphabet of size l, the construction of a balanced WT needs O a dlog l= ffiffiffiffiffiffiffiffiffi ffi log a p e À Á time (Munro, Nekrich & Vitter, 2016) and requires nH 0 (A) + o(a log l) bits when bit vectors are zero-order compressed .
The usual full-text indexes such as the FM-index (Ferragina & Manzini, 2000) or the LZ-index (Kärkkäinen & Ukkonen, 1996) do not index labeled texts. Some recent researches allow bidimensional range queries: Arroyuelo et al. (2015) represent an XML tree structure as a compressed bit vector and allow queries on structured patterns.
We focus here on non-overlapping labels and look for ways to efficiently store such labels along with sequences as well as to query them. The basic idea is that consecutive labels in A are often the same: we will thus store compressed label sequences.

TL-index: indexing labels over a text
Given a labeled text (T, A), we define the TL-index as, using a FM-index built on a BWT U to index the text, a bit vector B A marking the positions in the text where the labels change, and a WT W A indexing a compressed label sequence ( Fig. 2A).  Figure 1 The text T = AACAGC$ATCAAC$AGCTTT$, with three sequences, is labeled with the label Full-size  DOI: 10.7717/peerj-cs.148/ fig-1 BWT U. Let U = u 0 u 1 : : : u n -1 be the BWT of T = t 0 t 1 : : : t n -1 . As usually done, the FM-index samples every log 1+3 n values of a suffix array to retrieve the text positions of any occurrence . ε 001000010010 Wavelet Tree W A . Let A′ =〈a i j B A [i] = 1〉 . A′ = a′ 0 a′ 1 : : : a′ a -1 is called the compressed label sequence. It is a subsequence of A, of length a, containing only the changing labels according to the positions in A. The compressed label sequence A′ is stored using a WT W A . W A is defined as (Q, q 0 ), where Q is the set of nodes of W A and q 0 ∈ Q is the root node. Each node q ∈ Q is q = (q.val, q.label, q.left, q.right, q.parent), where q.label ∈ L ∪ {3} is the label of q and q.parent is the parent node of q (q 0 .parent being the null value ⊥). Both q.left and q.right are child nodes in Q ∪ {⊥}. A leaf is a node q where q.right and q.left are ⊥. Each leaf q is linked to a label q.label ∈ L ∪ {3}. The l + 1 leaves are exactly the labels of L ∪ {3} and we define leaf (q.label) = q. On a leaf q, we have q.val = ⊥. Let q be a non-leaf node: we have q.label = q.val is the bit vector rooted at q in the WT. We explain in "Shaping the WT for a Label Hierarchy" how W A can be further shaped depending on a label hierarchy.
W A is part of the index. This WT is used to answer efficiently bidimensional range queries where the labels and the label positions are the two dimensions . A balanced WT has a height of log l, with l, the number of leaves. The accessor W 〈i〉returns a′ i , in O(log l) time. This is a classical query within a WT. Given a label L x ∈ L, the function selectW 〈L x , i〉gives the position of the ith L x label in A′ in O(log l) time. The accessor W -1 〈L x 〉gives the list of positions in A′ where a′ i = L x . It runs in O(log l Â occ) time, with occ the number of occurrences of L x in A′.

TL BW -index: indexing labels in the order of the BWT
The BWT tends to store text repetitions consecutively. As those repetitions may have the same labels, it would be interesting that the labels benefit from the BWT reordering. Hence, labels can also be stored in the order of U.
Given a labeled text (T, A), the TL BW -index is defined as (U, B D , W D ) (Fig. 2B). The BWT U is built in the same way as the TL-index. Let D = d 0 d 1 : : : d n -1 the labels in the order of U. The bit vector B D of size n is such that B D [0] = 1, and, for i ! 1, The WT W D now indexes the compressed label sequence D′. The TL BW -index will be slower to build than the TL-index as it needs D. On the other side, as it is aware of the order of letters in the BWT, it will be able to support faster text/label queries.

Queries
The indexes allow the following classical queries.
label(t i ) (also called access(a i ))-Which label is on the letter t i ? This query is done is O(log l) time in the TL-index, and in O(log 1+3 n + log l) time in the TL BW -index since it has to convert positions in U order (see Fig. 3). findP(P) P-Which are the occurrences of a pattern P ? It is solved with the FM-index alone (Ferragina & Manzini, 2000). This query runs in O(jPj + occ Â log 1+3 n) time in both indexes, where occ is the number of occurrences of P in T.
The three previous queries are well known in text indexes. The two next queries search for a pattern and a label at the same time.
countPL(P, L x )-How many text positions are labeled L x and begin an occurrence of a pattern P ? As in the findL(L x ) query, the occurrences of P are found in U, in the positions from i to j.

TL-index:
We translate all these occ occurrences to have the corresponding positions in the text. For each of them, we run the query label(t i ). The total time is O(jPj + occ (log n 1+3 + log l)).
TL BW -index: See Algorithm 1. i and j correspond to the positions i′ = rank(1, i, B D ) and j′ = rank(1, j, B D ) in the root q 0 of W D . We then use an accessor customized from the rangeLocate function of , simulating a two-dimensional range query on This accessor first traverses the WT from the root to the leaf L x and then traverses it back to the root to find the positions in q 0 .val which match the label in the given range. For every position found, we find the corresponding positions in B D and expand them to the following 0-bits in B D . This query runs in O(jPj + jZj Â log l) time.
findPL(P, L x )-Which text positions are labeled L x and begin an occurrence of a pattern P?
TL-index: This query is exactly the same as countPL(P, L x ). if i′ > j ′ then return 0 As jZj y occ, the countPL() and the findPL() queries may thus be faster on the TL BW -index, jZj depending of the compression B D can do on Y.
Note that the countPL(P, L x ) query could be faster if the WT was directly built on the labels, without the intermediate bit vector B D (B A ): the answer would be known while reaching the leaf of the WT in O(jPj + log l) time. We chose to favor the execution time of the findPL(P, L x ) query, as well as the size of the structure (when A can be compressed).
The findPL(P, L x ) can vary to find the patterns P which have the label L x on the ith position of the pattern, or in any pattern's position. Adapting this queries is easy in the TL-index as we find the positions of the pattern in the BWT, translate them in the text order and then find the label of the ith position following each of them. In the TL BW -index, we find the patterns' positions in the BWT, access to the ith letter (we need to sample the BWT to read the patterns in the forward way), and find the label as usual.
To have the label in any pattern's position, in the TL-index we need to find a label L x between the first and last letter of the pattern (with only two access in the WT) but in the TL BW -index, we look for the label of all the pattern's letters.

Construction and space
We recall that the text (T, A) is of length n and is labeled with l unique labels. As defined above, the indexes store U in nH k (T) + o(n) bits. The TL-index stores the bit vector with rank and select capabilities in nH 0 (B A ) + o(n) bits. The size of W A depends on the compressed label sequence A′, of length a. W A takes a H 0 (A′) + o(a log l) bits. Similarly, the TL BW -index stores B D in nH 0 (B D ) + o(n) bits and W D takes d H 0 (D′) + o(d log l) bits, where d is the length of D′. The BWT can be built in linear time while using little space (Belazzougui et al., 2016;Munro, Navarro & Nekrich, 2017). B A is built while reading A in O(n) time. To make B D , we need to read the labels in the order of the original data file in O(n) time. To make W A , we find the occurrence of each label, corresponding to a 1-bit in B A , in O(a) time. Then we form the shape of W A in O(l). The labels corresponding to a 1-bit are extracted to make the WT's root q 0 . For each node containing at least two labels, we separate them by following the shape previously calculated, in O a log l= ffiffiffiffiffiffiffiffiffi ffi log a p AE Ç À Á . We build W D the same way. The TL-index has thus a size of nH k (T) + nH 0 (B A ) + a H 0 (A′) + o(n log l) bits, assuming s = O(l), and is built in O n þ l þ a log l= ffiffiffiffiffiffiffiffiffi ffi log a p AE Ç À Á time. The TL BW -index has a size of nH k (T) + d H 0 (D′) + o(n log l) bits and is built in O n þ l þ d log l= ffiffiffiffiffiffiffiffiffi ffi log d p AE Ç À Á time.

Shaping the WT for a label hierarchy
Labels may be further organized into a label hierarchy, given an additional set F = {F 0 , : : : , F f -1 } of label families (Fig. 5A). Both TL-and TL BW -indexes can be adapted: The WT W (either W A or W D ) will be shaped according to the hierarchy, and internal nodes q of W may have non-empty q.label values belonging to F. For example, in Fig. 5B, one can set on either index q 1 .label = L 1 , where L 1 is the label family gathering L 1.1 , L 1.2 and L 1.3 . The findL() and findPL() queries naturally extend to label families. With the hierarchy depicted in Fig. 5, findL(L 1 ) has to find the positions that have a label L 1.1 , L 1.2 or L 1.3 . Such a query does not need to be iterated on every label of the family L 1 , but rather directly starts at the corresponding internal node (q 1 on Fig. 5B).
Shaping W for a label hierarchy may increase the height w of the WT to O(l) in the worst case. To have a lower space consumption and a faster average query time, one can give a Huffman shape to W (Huffman, 1952). A leaf which corresponds to a frequently used label will be higher in the tree than a leaf of a rarely used label. Depending on the label hierarchy, the average height of the tree is H 0 (A′) in the best case while in the worst case it could be O(l). If no label hierarchy is given, the average height w will be H 0 (A′) (Mäkinen & Navarro, 2005).

HT-index: a baseline index
We compared the TL-and TL BW -indexes with a baseline solution called HT-index, indexing the text T with a BWT. The labels are stored in a map linking each label to the list of its (start,end) positions. We also store the labels in the text order with the compressed bit vector B A and, stored in plain form, A′. This enables the findL(L x ) query in O(y), where Y′ is the list of pairs (start,end) which represent the occurrences and y = jY′j. Note that Y (in the TL-index) and Y′ represent the same information as the labels are stored in the text order in both indexes. The label(t i ) query runs in O(1) time. This solution is not space-efficient with repeated labels: it needs nH k (T) + nH 0 (B A ) + a + l′ + o(n) bits, where a is the size of A′ and l′ the number of labeled factors. The query times are summarized in Table 1.

Evaluation procedure
The three indexes were implemented in C++. We used the SDSL-Lite library (Gog et al., 2014) to build the bit vectors and the WT. We used the RopeBWT2 library (Li, 2014), which builds a BWT in O(n log n) time on small DNA sequences, as it is very efficient for storing and appending sequences corresponding to our projected application (Cox et al., 2012). As RopeBWT2 does not sample the suffix array, we iterate over the text until we find a $ symbol. To have results close to the usual FM-index sampling in O(log 1+3 n) steps, we use sequences of length 50, which is similar to the sampling distance usually chosen in practice. The real files have longer sequences, thus longer sampling distances. The queries relying on the BWT will be slower and therefore cannot be compared between real and simulated files. We build the B A (or B D ) bit vectors, compress them using the rrr_vector class of SDSL-Lite, and finally build W A (or W D ) using a shape we added in SDSL-Lite, which depends of the label hierarchy. We evaluated the build time, the index size, as well as the run time of three of the queries detailed in "Queries"-findP(P) behaving similarly in all indexes and countPL(P, L x ) being very similar to findPL(P, L x ). The three indexes were tested on various datasets of labeled texts, each with 100 M characters (Table 2). Datasets and code are available at http://www.vidjil.org/data/#2018-peerjcs.
Simulated files with random sequences and random labels. All sequences and labels are random (d ∼ 0.8n). Simulated files, with random sequences but fixed labels. Here a given label is associated to the same pattern, possibly with some variation, and we alter the proportion of labeled letters (5-100%), the variation in the label's pattern (0-50%, more variations giving random labels), the number of unique labels (10-1,000), the length of the labels (5-100 letters). The dataset has 546 files, two of those files are shown in Table 2. Genomic sequences with immunologic labels. A person's immune system can be described with a set of labeled DNA sequences with V(D)J recombinations. The dataset, detailed below, uses 838 K labels from 355 unique labels, with d ∼ 0.26n.

A dataset of DNA sequences with immunologic labels
The adaptive immunity is the mechanism thanks to which our body defends itself against infections. When B and T-cells, or lymphocytes, are under maturation, their DNA undergo a recombination which allows several billions possibilities from a register of a thousand genes (Tonegawa, 1983). For example, the V(D)J recombination V4 Ã 02 4/ACGT/0 J1 Ã 12 means that the DNA sequence is composed from the V4 Ã 02 gene without the four last letters, then the ACGT sequence, then the J1 Ã 12 gene.
A person's immune system can thus be described with a set of labeled DNA sequences encoding V(D)J recombinations (Fig. 6). These sequences can be sampled by nextgeneration sequencing with bioinformatics analysis (Bystry et al., 2016, Duez et al., 2016. The tested DNA sequences come from patients 09, 12, 14 and 63 from a public dataset on a study on acute lymphoblastic leukemia (Salson et al., 2017). They have 100 M letters Table 1 Query time complexities for HT-index, TL-index and TL BW -index.

Requests
HT-index TL-index TL BW -index Note: Note that we have jZj y occ p . The label(i) and findL(L x ) queries are faster in the HT-index and the TL-index as the HT-index needs a sampling time. However, the countPL(P, L x ) and findPL(P, L x ) are faster in the HT-index. and 838 K labels from 355 unique labels, making a 117 MB file. Each DNA sequence has between 100 and 350 letters and two or three labels, each label having a size between 5 and 200 letters (Fig. 7). For a given label, the labeled subsequences may vary up to 15% due to amplification and sequencing errors.

Results
Index sizes. The Table 1 shows the results. As expected, the size of U, B (either B A or B D ) and W (either W A or W D ) grows linearly with the number of indexed elements (data not shown). The TL-index is the smallest index, and the TL BW -index is generally slightly larger. The compression is directly related to a and d. The file with random labels (d = 0.77n) is hard to compress, whereas the files with a low d/n ratio give a 2Â to 7Â compression. Figure 8 further shows how these sizes vary: As expected, the indexes are larger when there are less consecutive identical labels in T or in U, thus more 1-bits in B. Notes: All files have 100 M characters, and differ by the number of total and unique labels ("Lab (t/u)") and their size ("Lab. avg size"), by the ratio of labeled letters ("Lab. letters"), and by the variation between sequences labeled by the same label ("Variation"). Queries use patterns P with three letters. Times were averaged on 1 M launches (label()) or at least five launches (other queries). Times for findL(L) are reported per letter. Best or close-to-the-best results are in bold. Note that when there are more labeled letters, the text is more similar (as labels are placed on similar substrings), hence a decrease in the BWT size (A, C). W increases while the number of unique labels increases (D), the height of W increasing logarithmically with the number of unique labels. Build time. Most of the build time of TL BW -index is devoted to build D′. For TL-index, the building of U takes most of the total time. Queries. The label() query is faster on the HT-index. As expected, the TL BW -index needs more time for label() and findL(), as it needs to translate positions from the BWT back to the text. Note that locating the positions in the text takes about the same time as label(t i ) in the TL-index. However, for the complex findPL(P, L) query, the TL BW -index is the fastest solution because the position translation is only done on the letters which have both the label and the pattern. For the TL-index and HT-index, the actual time of the findPL(P, L) queries is more affected by the number of pattern occurrences than the number of final occurrences (between 0 and 100 K depending on the file).

A G G T C A A T A C G A T G A C T G G G G T C A G C T C A T A C G T C A G G A G G
On the genomic dataset, the sequences are longer: The TL BW -index suffers here even more from the missing suffix array sampling of the implementation for queries label() and findL(). However, on the findPL(P, L) query, the other indexes are penalized due to the sparsity of the sampling, bringing a more than 30Â difference with TL BW -index.

CONCLUSION
The TL-index and TL BW -index store a labeled text and allow efficient queries. They can be constructed from files of some MB in a few seconds. Experiments confirm that the indexes stay small when the text is redundant (thus a smaller U), when each label describes a pattern with few variations (many 0-bits in B, thus a smaller W), and when few letters are labeled (thus a small W). However, the TL-index and TL BW -index are robust even to nonredundant text with almost random labels everywhere. The TL BW -index needs more memory space than the TL-index but is more efficient in combined label/pattern queries. Those structures might be used on any labeled data, such as DNA sequences, but also on natural language texts or on music sheets with some semantic annotations.
Perspectives include improvement of the implementation, with label families queries or parameterizing the distance between samples in the FM-index to offer a space-time trade off. Within SDSL we could use the sd_vector bit vector instead of the rrr_vector bit vector which should improve space consumption when the bit vectors are very sparse. However, this would only minimally improve the global space consumption of the index. We plan to use one of the indexes in a clone database for hematological malignancies: It will allow comparison of V(D)J recombinations between different samples of a patient or between several patients.