Indexing labeled 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).

BWT U. Let U = u₀ u₁ … u_{n − 1} be the BWT of T = t₀ t₁ … t_{n − 1}. As usually done, the FM-index samples every log^1+ε n values of a suffix array to retrieve the text positions of any occurrence (Navarro & Mäkinen, 2007).

Bit vector B_A. Let B_A a compressed bit vector of length n such that B_A[0] = 1, and, for i ≥ 1, B_A[i] = 0 if a_i = a_{i − 1}, and otherwise B_A[i] = 1.

Wavelet Tree W_A. Let A′ = 〈a_i | B_A[i] = 1〉. A′ = a′₀ a′₁ … 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₀), where Q is the set of nodes of W_A and q₀ ∈ 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 ∪ {ε} is the label of q and q.parent is the parent node of q (q₀.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 ∪ {ε}. The l + 1 leaves are exactly the labels of L ∪ {ε} 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 (Mäkinen & Navarro, 2007). 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⁻¹〈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₀ d₁ … 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, B_D[i] = 0 if d_i = d_{i − 1}, and otherwise B_D[i] = 1. Let D′ = 〈d_i | B_D[i] = 1〉. D′ = d′₀, d′₁…d′_{d − 1} is a compressed label sequence of length d, subsequence of D. 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+ε 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(|P| + occ × log^1+ε n) time in both indexes, where occ is the number of occurrences of P in T.

• findL(L_x)—Which text positions are labeled L_x? The query runs in O(y × log l) time in the TL-index, and in O(y(log l + log^1+ε n)) time in the TL_BW-index, where Y = W⁻¹〈L_x〉 and y = |Y|. See Fig. 4.

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(|P| + occ (log n^1+ε + 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₀ of W_D. We then use an accessor customized from the rangeLocate function of (Mäkinen & Navarro, 2007), simulating a two-dimensional range query on [L_x, L_x] × [i′, j′] in W_D to give the set of positions Z = {z | a′_z = L_x and i′ ≤ z ≤ j′ } in O(|Z| × log l) time. 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₀.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(|P| + |Z| × 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).

  • –

    TL_BW-index: We use the countPL(P, L_x) query detailed in Algorithm 1, replacing the counter cnt with a list Y holding the matching positions in U. The positions are converted at the end in the text order. This query runs in O(|P| + |Z| × log l + y × log^1+ε n) time, with |Y| = y.

As |Z| ≤ y ≤ occ, the countPL() and the findPL() queries may thus be faster on the TL_BW-index, |Z| 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(|P| + 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₀(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₀(A′) + o(a log l) bits. Similarly, the TL_BW-index stores B_D in nH₀(B_D) + o(n) bits and W_D takes d H₀(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₀. For each node containing at least two labels, we separate them by following the shape previously calculated, in $O\left(a\lceil \log l/\sqrt{\log a}\rceil \right)$. We build W_D the same way. The TL-index has thus a size of nH_k(T) + nH₀(B_A) + a H₀(A′) + o(n log l) bits, assuming σ = O(l), and is built in $O\left(n+l+a\lceil \log l/\sqrt{\log a}\rceil \right)$ time. The TL_BW-index has a size of nH_k(T) + d H₀(D′) + o(n log l) bits and is built in $O\left(n+l+d\lceil \log l/\sqrt{\log d}\rceil \right)$ time.

Shaping the WT for a label hierarchy

Labels may be further organized into a label hierarchy, given an additional set F = {F₀, …, 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₁.label = L₁, where L₁ 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₁) 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₁, but rather directly starts at the corresponding internal node (q₁ 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₀(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₀(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 = |Y′|. 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₀(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+ε 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 next-generation 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 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. 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).

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.