Fast and Memory-Efficient BWT Construction of Repetitive Texts Using Lyndon Grammars

The Burrows-Wheeler Transform ($\mathrm{𝖡𝖶𝖳}$) [13] is one of the most important data structures in sequence analysis and is applied in a wide range of fields from data compression to bioinformatics. It is obtained by assigning the $i$th symbol of the $\mathrm{𝖡𝖶𝖳}$ to the last character of the $i$th lexicographically smallest conjugate of the input text. The $\mathrm{𝖡𝖶𝖳}$ can be computed in linear time and space. The run-length compressed $\mathrm{𝖡𝖶𝖳}$ ($\mathrm{𝖱𝖫𝖡𝖶𝖳}$) uses the fact that the number of equal-character runs in the $\mathrm{𝖡𝖶𝖳}$ is small compared to the text length when the text is repetitive [13]. Mantaci et al. [32] extended the notion of the $\mathrm{𝖡𝖶𝖳}$ to string collections: the extended $\mathrm{𝖡𝖶𝖳}$ ($\mathrm{𝖾𝖡𝖶𝖳}$) of a string collection consists of the last characters of the strings in $ℳ$ arranged in infinite periodic order (see Section 2 for definitions).

The $\mathrm{𝖡𝖶𝖳}$ serves as the basis of several self-indexes, i.e., data structures supporting access to the text as well as pattern matching queries. One of the most successful such indexes is the FM-index [20], which is used e.g. in the important bioinformatics programs BWA [30] and Bowtie [28]. However, the FM-index uses space proportional to the size of the dataset and is thus not applicable for huge datasets that vastly exceed the size of main memory. Such datasets often contain data from thousands or millions of individuals from the same species and are therefore extremely repetitive. The r-index [22] is based on the $\mathrm{𝖱𝖫𝖡𝖶𝖳}$ and is the first text index of size proportional to the number $r$ of runs in the $\mathrm{𝖡𝖶𝖳}$. It can be constructed in space proportional to $r$ from the $\mathrm{𝖱𝖫𝖡𝖶𝖳}$, which raises the problem of constructing the $\mathrm{𝖱𝖫𝖡𝖶𝖳}$ of huge datasets with small working memory. For this reason, various tools have emerged that exploit the repetitiveness of these datasets.

The Lyndon grammar of a text is a straight-line program (SLP), i.e., a context-free grammar in Chomsky normal form that generates a single text, where each symbol corresponds to a node in the Lyndon tree. It was introduced in [25] as the basis for a self-index. In this paper, we describe a method for efficiently computing several $\mathrm{𝖡𝖶𝖳}$ variants of a text or a text collection via the input’s Lyndon grammar and show that our method surpasses other tools with respect to time or memory used, especially when using multiple threads.

For $i,j\in {\mathbb{N}}_0$ we denote the set $\{k\in {\mathbb{N}}_0:i\le k\le j\}$ by the interval notations $[i..j]=[i..j+1)=(i-1..j]=(i-1..j+1)$. For an array $𝖠$ we analogously denote the subarray from $i$ to $j$ by $𝖠[i..j]=𝖠[i..j+1)=𝖠(i-1..j]=𝖠(i-1..j+1)=𝖠[i]𝖠[i+1]\mathrm{\dots }𝖠[j]$. We use zero-based indexing, i.e., the first entry of the array $𝖠$ is $𝖠[0]$.

A string $S$ of length $n$ over an alphabet $\mathrm{\Sigma }$ is a sequence of $n$ characters from $\mathrm{\Sigma }$. We denote the length $n$ of $S$ by $|S|$ and the $i$th symbol of $S$ by $S[i-1]$, i.e., strings are zero-indexed. In this paper we assume any string $S$ of length $n$ to be over a totally ordered and linearly sortable alphabet (i.e., the characters in $S$ can be sorted in $𝒪(n)$). Analogous to arrays we denote the substring from $i$ to $j$ by $S[i..j]=S[i..j+1)=S(i-1..j]=S(i-1..j+1)=S[i]S[i+1]\mathrm{\dots }S[j]$. For $j>i$ we let $S[i..j]$ be the empty string $\epsilon$. For two strings $u$ and $v$ and an integer $k\ge 0$ we let $uv$ be the concatenation of $u$ and $v$ and denote the $k$-times concatenation of $u$ by $u^k$. A string $S$ is primitive if it is non-periodic, i.e., $S=w^k$ implies $w=S$ and $k=1$. The suffix $i$ of a string $S$ of length $n$ is the substring $S[i..n)$ and is denoted by ${\mathrm{suf}}_i(S)$. Similarly, the substring $S[0..i]$ is a prefix of $S$. A suffix (prefix) is proper if $i>0$ (). For two possibly empty strings $u$ and $v$, $uv$ is a conjugate of $vu$.

We assume totally ordered alphabets. This induces a total order on strings. Specifically, we say a string $S$ of length $n$ is lexicographically smaller than another string $T$ of length $m$ if and only if there is some $\mathrm{\ell }\le \min \{n,m\}$ such that $S[0..\mathrm{\ell })=T[0..\mathrm{\ell })$ and either or and , and write in this case. A non-empty string $S$ is in its canonical form if and only if it is lexicographically minimal among its conjugates. If $S$ is additionally strictly lexicographically smaller than all of its other conjugates, $S$ is a Lyndon word. Equivalently, $S$ is a Lyndon word if and only if $S$ is lexicographically smaller than all its proper suffixes [18].

In the following, we use abbabcbcabb as our running example.

The Lyndon factorization of our running example is $\text{abbabcbc},\text{abb}$.

The standard factorization of the first Lyndon factor of our running example is $(\text{abb},\text{abcbc})$ as abcbc is Lyndon and no longer proper suffix is Lyndon.

Figure 1 illustrates the following definitions.

For instance, and $\text{abb}{>}_{\omega }\text{ab}{>}_{\omega }\text{aba}$ since $\text{abb}\mathrm{\cdots }{>}_{\mathrm{𝑙𝑒𝑥}}\text{abab}\mathrm{\cdots }{>}_{\mathrm{𝑙𝑒𝑥}}\text{abaab}\mathrm{\dots }$.

Note that, by definition, the $\mathrm{𝖡𝖡𝖶𝖳}$ of a string $S$ is the $\mathrm{𝖾𝖡𝖶𝖳}$ of the Lyndon factors of $S$. Similarly, the $\mathrm{𝖾𝖡𝖶𝖳}$ of a multiset $ℳ$ of strings is the $\mathrm{𝖡𝖡𝖶𝖳}$ of the concatenation of the canonical forms of the strings in $ℳ$ arranged in lexicographically decreasing order [3] (the Lyndon factors of the resulting string are conjugates of the roots of the strings in $ℳ$).¹¹1The $\mathrm{𝖾𝖡𝖶𝖳}$ was originally defined for sets of primitive strings [32]. This limitation is unnecessary [10, 40].

For a primitive string $S$, its $\mathrm{𝖡𝖶𝖳}$, the $\mathrm{𝖾𝖡𝖶𝖳}$ of $\{S\}$, and the $\mathrm{𝖡𝖡𝖶𝖳}$ of the smallest conjugate of $S$ are identical. Consequently, the $-$\mathrm{𝖡𝖶𝖳}$ of a string $S$ (i.e., $\mathrm{𝖡𝖶𝖳}(S\$)$, where $\$$ is smaller than all characters in $S$) commonly computed via the suffix array of $S$, is identical to $\mathrm{𝖡𝖶𝖳}(\$S)=\mathrm{𝖡𝖡𝖶𝖳}(\$S)=\mathrm{𝖾𝖡𝖶𝖳}(\{\$S\})$.

Consider a Lyndon grammar $𝒢$ where the symbols are sorted lexicographically, i.e., . In this section we examine such Lyndon SLPs and find properties that enable sorting any Lyndon SLP lexicographically in linear time (Section 4) and give rise to an online algorithm for constructing Lyndon SLPs (Section 6.3). Note that the grammar shown in Figure 2 is lexicographically sorted.

For example, we have ${𝒫}_L(X_6)=\{X_6,X_5,X_1\}$, ${𝒫}_R(X_6)=\{X_6,X_8,X_9\}$ and ${𝒫}_L^{-1}(X_1)=\{X_1,X_2,X_3,X_4,X_5,X_6\}$ for the grammar shown in Figure 2.

Let $r(X_i)$ be the largest symbol whose derivation tree has $X_i$ on the leftmost path, i.e., $r(X_i)=\max \left({𝒫}_L^{-1}(X_i\right))$. Similarly, let $l(X_i)$ be the smallest symbol whose derivation tree has $X_i$ on the leftmost path, i.e. $l(X_i)=\min \left({𝒫}_L^{-1}(X_i\right))$. For example, $l(X_1)=X{}_1$ and $r(X_1)=X{}_6$ in the grammar from Figure 2. In the remainder of this section, we show that all symbols $X_j$ in the interval $[l(X_i)..r(X_i)]$ generate strings with $[X_i]$ as prefix and even satisfy $X_j\in {𝒫}_L(X_i)$. This then implies that these intervals form a tree-structure.

Given a Lyndon SLP $𝒢$ with roots $r_1,\mathrm{\dots },r_k$, we want to rename the symbols such that holds if and only if . This is possible in $𝒪\left(|𝒢\right|)$ using a similar principle to what is used for computing the Lyndon grouping in [2, 40]: consider a symbol $X_i$. By Lemma 17 and simple properties of Lyndon words, the symbols in $𝒟(X_i)\setminus {𝒫}_L(X_i)$ are lexicographically greater than $X_i$. Therefore, when considering the symbols in lexicographically decreasing order we can “induce” the position of a symbol $X_i\to X_a{X}_b$ upon encountering $X_b$: the symbols in ${𝒫}_L^{-1}(X_i)$ must be the lexicographically largest symbols with $X_a$ as prefix symbol that have not yet been induced.

Algorithm 1 shows the procedure. Throughout, $𝖠[i]$ is either the $i$th symbol in the lexicographical order or $\perp$, and for each inserted $X_i$, $𝖦[i]$ is the index in $𝖠$ of the smallest child of $X_i$ that has been inserted into $𝖠$ (or $r(X_i)+1$ if no child has been inserted yet). In the first for-loop, all terminal symbols are inserted into $𝖠$, and in the second for-loop, the remaining symbols are induced. Consequently, after the second for-loop, $𝖠$ contains the lexicographic order of the symbols, and $𝖦$ is (almost) it’s inverse, i.e. $𝖦[𝖠[i]]=i+1$.

There are two operations in Algorithm 1 that are not immediately obvious, namely finding $\left|{𝒫}_L^{-1}(X_i\right)|$ and iterating over all $X_j$ with $X_j\to X_a{X}_{𝖠[i]}$ for a given $𝖠[i]$. For the first, note that the values $\left|{𝒫}_L^{-1}(X_i\right)|$ can be trivially computed in linear-time using the recurrence $\left|{𝒫}_L^{-1}(X_a\right)|=1+{\sum }_{X_i\to X_a{X}_b}\left|{𝒫}_L^{-1}(X_i\right)|$. Secondly, iterating efficiently over all $X_j$ with $X_j\to X_a{X}_{𝖠[i]}$ can be achieved using a linear-time preprocessing step where we collect for each symbol $X_i$ the symbols which have $X_i$ as second symbol on the right-hand side.

The following theorem follows. Its proof is omitted due to space constraints.

In this section, we describe an algorithm that derives the $\mathrm{𝖡𝖡𝖶𝖳}$ from the Lyndon SLP. Specifically, combining the second phase of GSACA [2, 40] with run-length encoding on the lexicographically sorted Lyndon SLP results in a very efficient algorithm on real data.

The following lemmas and corollary establish a relationship between the lexicographical order of suffixes and the lexicographical order of certain symbols of the Lyndon grammar.

As shown in [1], the simple folklore algorithm for constructing the Lyndon array can be used to construct the Lyndon forest. With a fitting lookup data structure, we can trivially construct the Lyndon SLP instead. We essentially use Algorithm 4.

There are two steps that have to be explained: First, one must be able to compare $[X_t]$ and $[X_c]$ lexicographically, and secondly, one must be able to find $X_{c^{\prime }}$ with $X_{c^{\prime }}\to X_c{X}_t$ (this is also called the naming function). For the former, we store for each symbol on the stack a fixed-length prefix of the generated string. In many cases, this is sufficient for determining the lexicographical order. Possibilities to handle the other cases are explained in Sections 6.1, 6.2 and 6.3. The following paragraphs will deal with the naming function.

For the naming function, we use a hash table to find the symbol names. This provides (expected) constant time lookup/insertion.⁶⁶6Note that this is also possible in $𝒪(\log \log g)$ deterministic time with at most $2g+g\log g+o(g\log g))$ bits of memory, where $g$ is the number of symbols of the grammar [42]. In practice however, $N$ hash table lookups are unnecessarily slow, especially because we assume the input to be repetitive. In particular, having fewer but longer keys is cache friendlier and thus faster on modern computer architectures.

Specifically, we decrease the number of hash table accesses by partitioning the nodes of the Lyndon forest into heavy and light nodes and querying a second hash table $H$ for each heavy node. The naming function is then only used when such a query to $H$ does not return an answer. Leaves of the Lyndon forest are considered heavy. An inner node $v$ is considered heavy, if and only if its number of immediate heavy descendants (heavy descendants for which only light nodes occur on the path to $v$) exceeds a fixed constant $n_{\mathrm{𝑡ℎ𝑟𝑒𝑠}}>1$ (see Figure 4).⁷⁷7In our experiments, $n_{\mathrm{𝑡ℎ𝑟𝑒𝑠}}=31$.

For a heavy node $v$ of the Lyndon forest we query $H$ with the sequence of immediate heavy descendants $v_1,\mathrm{\dots },v_k$ of $v$. This means that we only store heavy nodes on the stack, for a light node we instead store all immediate heavy descendants.

Different occurrences of the same symbol represent the exact same Lyndon tree and therefore correspond to the exact same sequence of immediate heavy descendants. Consequently, if a subtree isomorphic to the subtree rooted at a heavy node has occurred previously, the query to $H$ is sufficient to resolve it. Otherwise, we use the naming function to find the correct symbol and insert the sequence of immediate heavy descendants with it into $H$.

Note that each heavy node has at most $2n_{\mathrm{𝑡ℎ𝑟𝑒𝑠}}$ immediate heavy descendants because it has exactly two children (or zero), each of which is either light (and thus has at most $n_{\mathrm{𝑡ℎ𝑟𝑒𝑠}}$ immediate heavy descendants) or heavy itself. This implies that the total memory usage increases by a factor of at most $𝒪(n_{\mathrm{𝑡ℎ𝑟𝑒𝑠}})$.

The source code of our implementation is publicly available.⁸⁸8https://gitlab.com/qwerzuiop/lyndongrammar

We compare our $\mathrm{𝖡𝖶𝖳}$ algorithms for single texts with the programs Big-BWT ⁹⁹9https://gitlab.com/manzai/Big-BWT, last accessed: 22.04.2025, git hash 944cb27 [11] and r-pfbwt ¹⁰¹⁰10https://github.com/marco-oliva/r-pfbwt, last accessed: 22.04.2025, git hash 1fea5c3 [41], as well as libsais.¹¹¹¹11https://github.com/IlyaGrebnov/libsais, last accessed: 22.04.2025, git hash a138159 The latter uses a modified version of the Suffix-Array Induced Sorting (SAIS) algorithm [37] to compute the $\mathrm{𝖡𝖶𝖳}$ and, since it is based on the currently fastest $\mathrm{𝖲𝖠}$ construction implementation for general real-world data [40], can be viewed as a lower bound for algorithms using the suffix array to compute the $\mathrm{𝖡𝖶𝖳}$.

For the $\mathrm{𝖡𝖶𝖳}$ of text collections, we compare with PFP-eBWT ¹²¹²12https://github.com/davidecenzato/PFP-eBWT, last accessed: 22.04.2025, git hash 4ca75ce [10], r-pfbwt [41], ropebwt3 ¹³¹³13https://github.com/lh3/ropebwt3, last accessed: 22.04.2025, git hash 36a6411 [29], grlBWT ¹⁴¹⁴14https://github.com/ddiazdom/grlBWT, last accessed: 22.04.2025, git hash f09e7fa [17] and CMS-BWT ¹⁵¹⁵15https://github.com/fmasillo/cms-bwt, last accessed: 22.04.2025, git hash 1099d07. Note that the speed and memory usage of CMS-BWT has improved massively since its publication in [35]. [35] (for CMS-BWT we used the first sequence in the collection as reference). All tool except the last one support multi-threading. Note that not all of these tools compute the same $\mathrm{𝖡𝖶𝖳}$ variant [14]. Also note that all algorithms based on PFP as well as grlBWT are semi-external, i.e., write/read some temporary data to/from disk.

As test data, we use up to 1000 human Chromosome 19 haplotypes from [12] ($\approx 6\cdot {10}^{10}$bp) and ${10}^6$ SARS-CoV-2 sequences ($\approx 3\cdot {10}^{10}$bp).¹⁶¹⁶16Downloaded from https://www.covid19dataportal.org on 28.08.2024. All experiments were conducted on a Linux-6.8.0 machine with an Intel Xeon Gold 6338 CPU and 512 GB of RAM. All programs were compiled with GCC 13.3.0. Before each test, the test file was scanned once to ensure it is cached by the kernel. The results can be seen in Figures 5, 6 and 7. The subscripts for the PFP-based algorithms indicate the used modulus.

For our programs, computing the $\mathrm{𝖽𝗈𝗅𝖤𝖡𝖶𝖳}$ using the suffix comparison method from Section 6.2 is generally the fastest. Computing the original $\mathrm{𝖾𝖡𝖶𝖳}$ is slightly slower because we need to find the smallest rotation of each input string. In the single-threaded case, for both the 1000 Chromosome 19 haplotypes and ${10}^6$ SARS-CoV-2 sequences, over $98\%$ of the time was spent constructing the grammar (for our fastest algorithm). Using the linear-time algorithm for constructing the Lyndon array from [7] to ensure expected linear time complexity slows our programs down by up to $60\%$. As expected, the suffix comparison method from Section 6.3 is much slower than our other methods. The increase in memory consumption of our program in the multithreaded case is due to the use of thread-safe hash tables and multiple sequences and stacks residing in main memory.

For the Chromosome 19 collections and a single thread, CMS-BWT is the fastest program (at the cost of high memory usage), followed by r-pfbwt (for larger cases) and our algorithms. For more threads, our program is always the fastest. Regarding the memory consumption, grlBWT uses the least amount of main memory for single-threaded processing, followed by our programs. For the SARS-CoV-2 sequences, our program is the fastest, especially with multiple threads, and for more than one thread also the most memory efficient.

We described an algorithm to compute the $\mathrm{𝖡𝖡𝖶𝖳}$ – and by extension the common $-$\mathrm{𝖡𝖶𝖳}$ and various versions of the $\mathrm{𝖾𝖡𝖶𝖳}$ – from the lexicographically sorted Lyndon grammar of a text or text collection. Furthermore, we gave an algorithm that lexicographically sorts a Lyndon grammar and discussed approaches to efficiently compute the Lyndon grammar of a text or text collection. We implemented the described algorithms and found that they outperform other current state-of-the art programs in terms of time or memory consumed (often both).