BWT construction and search at the terabase scale

Abstract

Motivation

Burrows–Wheeler Transform (BWT) is a common component in full-text indices. Initially developed for data compression, it is particularly powerful for encoding redundant sequences such as pangenome data. However, BWT construction is resource intensive and hard to be parallelized, and many methods for querying large full-text indices only report exact matches or their simple extensions. These limitations have hampered the biological applications of full-text indices.

Results

We developed ropebwt3 for efficient BWT construction and query. Ropebwt3 indexed 320 assembled human genomes in 65 h and indexed 7.3 terabases of commonly studied bacterial assemblies in 26 days. This was achieved using up to 170 gigabytes of memory at the peak without working disk space. Ropebwt3 can find maximal exact matches and inexact alignments under affine-gap penalties, and can retrieve similar local haplotypes matching a query sequence. It demonstrates the feasibility of full-text indexing at the terabase scale.

Availability and implementation

https://github.com/lh3/ropebwt3.

1 Introduction

Although millions of genomes have been sequenced, the majority of them were sequenced from a small number of species such as human, Escherichia coli and Mycobacterium tuberculosis. As a result, existing genome sequences are highly redundant. This is how Hunt et al. (2024) compressed 7.86 terabases (Tb) of bacterial assemblies, also known as AllTheBacteria, into 78.5 gigabytes (GB) after grouping phylogenetically related genomes (Břinda et al. 2024). The resultant compressed files losslessly keep all the sequences but are not directly searchable. Indexing is necessary to enable fast sequence search.

K-mer data structures are a popular choice for sequence indexing (Marchet et al. 2021). They can be classified into three categories. The first category does not associate k-mers with their positions in the database sequences. These data structures support membership query or pseudoalignment (Bray et al. 2016), but cannot reconstruct input sequences or report base alignment. Sequence search at petabase scale use all such methods (Edgar et al. 2022, Karasikov et al. 2024, Shiryev and Agarwala 2024). The second category associates a subset of k-mers with their positions. Upon finding k-mer matches, algorithms in this category go back to the database sequences and perform base alignment. Most aligners work this way. However, because the database sequences are not compressed well, these algorithms may require large space to store them. The last category keeps all k-mers and their positions. Algorithms in this category can reconstruct all the database sequences without explicitly storing them. Nonetheless, although positions of k-mers can be compressed efficiently (Karasikov et al. 2020), they still take large space. The largest lossless k-mer index consists of a few terabases (Karasikov et al. 2024).

Compressed full-text indices, such as FM-index (Ferragina and Manzini 2000), r-index (Gagie et al. 2018, Bannai et al. 2020, Gagie et al. 2020), and move index (Nishimoto and Tabei 2021), provide an alternative way for fast sequence search (Navarro 2022). The core component of these data structures is often Burrows–Wheeler Transform (BWT; Burrows and Wheeler 1994) which is a lossless transformation of strings. The BWT of a highly redundant string tends to group symbols in the original string into long runs and can thus be well compressed. When we supplement BWT with a data structure to efficiently compute the rank of a symbol in BWT, we can in theory count the occurrences of a query string in time linear to its length. FM-index further adds a sampled suffix array to locate substrings, while r-index uses an alternative method that is more efficient for redundant strings. Both of them support compression and sequence search at the same time.

BWT-based indices have been used for read alignment (Langmead et al. 2009, Li and Durbin 2009, Li et al. 2009), de novo sequence assembly (Simpson and Durbin 2012), metagenome profiling (Kim et al. 2016) and data compression (Cox et al. 2012). They have also emerged as competent data structures for pangenome data. Existing pangenome-focused tools (Ahmed et al. 2021, Rossi et al. 2022, Zakeri et al. 2024) use prefix-free parsing for BWT construction (Boucher et al. 2019). They require more memory than input sequences and are impractical for huge datasets. Although ropebwt2 developed by us can construct BWT in memory proportional to its compressed size and is fast for short strings (Li 2014), it is inefficient for chromosome-long sequences. OnlineRlbwt (Bannai et al. 2020) has a similar limitation. grlBWT (Díaz-Domínguez and Navarro 2023) is likely the best algorithm for constructing the BWT of similar genomes. It reduces the peak memory at the cost of large working disk space. In its current form, the algorithm does not support update to BWT. We would need to reconstruct the BWT from scratch when new genomes arrive. BWT construction for highly redundant sequences remains an active research area.

With BWT-based data structures, it is trivial to test the presence of string P in the index, but finding substring matches within P needs more thought. Learning from bidirectional BWT (Lam et al. 2009), we found a BWT constructed from both strands of DNA sequences supports the extension of exact matches in both directions (Li 2012). This gave us an algorithm to compute maximal exact matches (MEMs), which was later improved by Gagie (2024) for long MEMs. Bannai et al. (2020) proposed a distinct algorithm to compute MEMs without requiring both strands. Matching statistics (Chang and Lawler 1994) and pseudo-matching length (PML; Ahmed et al. 2021) have also been considered. All these algorithms find exact local matches only.

It is also possible to identify inexact matches. Kucherov et al. (2016) invented search schemes to find semi-global alignment. With the BWT-SW algorithm, Lam et al. (2008) simulated a suffix trie, a tree data structure, with BWT and performed sequence-to-trie alignment to find all local matches between a query string and the BWT of a large genome. We went a step further by representing the query string with its direct acyclic word graph (DAWG; Blumer et al. 1983) and performed DAWG-to-trie alignment. This is the BWA-SW algorithm (Li and Durbin 2010). Nonetheless, we later realized the formulation of BWA-SW had theoretical flaws and its implementation was closer to DAWG-to-DAWG alignment than to DAWG-to-tree alignment.

In this article, we will explain our solution to BWT construction and to sequence search at the terabase scale. Our contribution includes: (i) a reinterpreted BWA-SW algorithm that fixes issues in our earlier work (Li and Durbin 2010); (ii) its application in finding similar haplotypes; (iii) an incremental in-memory BWT construction implementation that scales to large pangenome datasets; (iv) a faster algorithm to find alignment positions with a standard FM-index of highly redundant strings.

2 Materials and methods

In a nutshell, ropebwt3 computes the partial multi-string BWT of a subset of sequences with libsais (https://github.com/IlyaGrebnov/libsais) and merges the partial BWT into the existing BWT run-length encoded as a B+-tree (Li 2014). It repeats this procedure until all input sequences are processed. The BWT by default includes input sequences on both strands. This enables forward-backward search (Li 2012) required by accelerated long MEM finding (Gagie 2024). Ropebwt3 also reports local alignment with affine-gap penalty using a revised BWA-SW algorithm (Li and Durbin 2010).

Ropebwt3 combines and adapts existing algorithms and data structures. Nonetheless, the notations here differ from our early work (e.g. from 1-based to 0-based coordinates) and from other publications. We will describe our methods in full for completeness.

2.1 Basic concepts

Let Σ be an alphabet of symbols (Table 1). Given a string P over Σ, $|P|$ is its length and $P[i]∈Σ, 0≤i<|P|$⁠, is the ith symbol in P, and $P[i,j)$ is a substring of length j—i starting at i. Operator “$°$” concatenates two strings or between strings and symbols. It may be omitted if concatenation is apparent from the context.

Suppose $T=(P0,P1,…,Pm−1)$ is an ordered list of m strings over Σ. $T≜P0$0P1$1⋯Pm−1$m−1$ is the concatenation of the strings in $T$ with sentinels ordered by $$0<$1<⋯<$m−1$⁠. For other ways to define string concatenation on ordered string lists or unordered string sets, see Cenzato and Lipták (2024).

For convenience, let $n≜|T|$ and $T[−1]=T[n−1]$⁠. The suffix array of T is an integer array S such that S(i), $0≤i<n$⁠, is the start position of the ith smallest suffix among all suffixes of T (Fig. 1a). The Burrows–Wheeler Transform (BWT) of T is a string B computed by $B[i]=T[S(i)−1]$⁠. All sentinels in B are represented by the same symbol “$$$” and are not distinguished from each other. $Σ′≜Σ∪{$}$ denotes the alphabet including the sentinel.

For $a∈Σ′$⁠, let $CB(a)≜|{i:B[i]<a}|$ be the number of symbols smaller than a and $rankB(a,k)≜|{i<k:B[i]=a}|$ be the number of a before offset k in B. We may omit subscription B when we are describing one string only. The last-to-first mapping (LF mapping) π is defined by $π(i)≜S−1(S(i)−1)$⁠, where $S−1$ is the inverse function of suffix array S. It can be calculated as $π(i)=πB(B[i],i)$⁠, where $πB(a,i)≜C(B[i])+rank(B[i],i)$⁠. As $B[π(i)]$ immediately proceeds $B[i]$ on T, we can use π to decode the ith sequence in B.

2.2 BWT construction

Libsais is an efficient library for computing the suffix array of a single string. It does not directly support a list of strings. Nonetheless, we note that T is a string over alphabet $Σ″={$0,$1,…,$m−1,A,C,G,T,N}$ with lexicographical order $$0<⋯<$m−1<A<C<G<T<N$⁠. We can use m + 5 non-negative integers to encode T and apply libsais. The suffix array derived this way will be identical to the suffix array of T.

For m that can be represented by a 32-bit integer and n represented by a 64-bit integer, libsais will need at least 12n bytes to construct the suffix array. It is impractical for n more than tens of billions. To construct BWT for large n, ropebwt3 uses libsais to build the BWT for a batch (up to 7 Gb by default) and merges it to the BWT of already processed batches (Algorithm 1). The basic idea behind the algorithm is well known (Ferragina et al. 2010) but implementations vary (Sirén 2016, Oliva et al. 2021). In ropebwt3, we encode BWT with a B+-tree (Fig. 2a). This yields $O(log r)$ rank query (line 13) and insertion (line 17), where r is the number of runs in the merged BWT. The bottleneck of the algorithm lies in rank calculation (line 8), which can be parallelized if $m2>1$⁠.

This online BWT construction algorithm does not use temporary disk space. The overall time complexity is $O(n log r)$⁠. The BWT takes $O(r log r)$ in space. The memory required for partial BWT construction with libsais depends on the batch size and the longest string. B+-tree is dynamic. Ropebwt3 optionally converts B+-tree to the fermi binary format (Fig. 2b; Li 2012) which is static but is faster to query and can be memory-mapped.

Ropebwt2 (Li 2014) uses the same B+-tree to encode BWT and has the same time complexity. However, it inserts sequences, not partial BWTs, into existing BWT. It cannot be efficiently parallelized for long strings. Note that independent of our earlier work, Ohno et al. (2018) also used a B+-tree to encode BWT. Its implementation (Bannai et al. 2020) is several times slower than ropebwt2 on 152 bacterial genomes from Li et al. (2024), possibly because ropebwt2 is optimized for the small DNA alphabet.

2.3 Suffix array interval and backward search

For a string $P∈Σ*$ (ie not including sentinels), let $occ(P)$ be the number of occurrences of P in T. Define $lo(P)$ to be the number of suffixes that are lexicographically smaller than P and $hi(P)≜lo(P)+occ(P)$⁠. $[lo(P),hi(P))$ is called the suffix array interval of P, or SA interval in brief. An SA interval is important if there exists P such that it is the SA interval of P. The SA interval of the empty string is $[0,n)$ where $n=|T|$⁠.

If we know the SA interval of P, we can calculate the SA interval of aP with: $lo(aP)=πB(a,lo(P))$ and $hi(aP)=πB(a,hi(P))$⁠. To calculate $occ(P)$⁠, we start with $[0,n)$ and repeatedly apply the equation above from the last symbol in P to the first. This procedure is called backward search.

2.4 Locating with FM-index

By definition of BWT, if $i∈[lo(P),hi(P)), P=T[S(i),S(i)+|P|)$⁠. We can thus locate an occurrence of P in the original string T. However, suffix array S takes $O(n log n)$ in space and may be too large to store explicitly. With an FM-index, we only store S(i) if and only if i is a multiple of s where s is a positive integer controlling the sample rate. We calculate the rest of S(i) using LF-mapping $π(·)$⁠.

A complication with multi-string BWT is that $π(i)$ does not point to the preceeding symbol when $B[i]=$$⁠. This is because sentinels are ordered during suffix sorting but are not distinguished from each other in BWT. We will have to store the rank of each string in T (array R in Algorithm 2). For convenience, we also keep the index of the string and the offset on the string instead of the offset on the concatenated string T.

Function Locate1 in Algorithm 2 provides a faster way to find one position in SA interval $[l,h)$⁠. A key observation is that if the interval contains a sentinel or there exists k such that $l≤ks<h$⁠, we can immediately locate one occurrence; if not, we can apply backward search repeatedly until an SA interval $[l′,h′)$ brackets a stored suffix array value. If T consists of $m′$ identical strings, we apply $s/min(h−l,m′)$ rounds of backward searches on average, usually much faster than the naive algorithm. Ropebwt3 implements a generalized Locate1 function that finds multiple occurrences.

2.5 Double-strand BWT

The definitions above are applicable to generic strings. With one BWT, we can only achieve backward search; forward search additionally requires the BWT of the reverse strings (Lam et al. 2009). Nonetheless, due to the strand symmetry of DNA strings, it is possible to achieve both forward and backward search with one BWT provided that the BWT contains both strands of DNA strings (Li 2012).

Formally, a DNA alphabet is $Σ={A,C,G,T,N}$⁠. $a¯$ denotes the Watson-Crick complement of symbol $a∈Σ$⁠. The complement of $$, A, C, G, T$⁠, and $N$ are $$, T, G, C, A$⁠, and $N$⁠, respectively.

For string P, $P¯$ is its reverse complement string. The double-strand concatenation of a DNA string list $T=(P0,P1,…,Pm−1)$ is $T˜=P0$0P¯0$1P1$2P¯1$3⋯Pm−1$2m−2P¯m−1$2m−1$⁠. The double-strand BWT (DS-BWT) of $T$ is the BWT of $T˜$⁠. We note that if P is a substring of $T˜$⁠, $P¯$ must be a substring and $occ(P)=occ(P¯)$⁠. The suffix array bidirectional interval (SA bi-interval) of P is a 3-tuple defined as $(lo(P),lo(P¯),occ(P))$⁠.

An SA bi-interval can be extended in both backward and forward directions. To calculate the SA bi-interval of aP, we can use the standard backward search to compute $lo(aP)$ and $occ(aP)$ from $lo(P)$ and $occ(P)$⁠. As to $lo(aP¯)$⁠, we note that $[lo(aP¯),hi(aP¯))⊂[lo(P¯),hi(P¯))$ because $aP¯=P¯°a¯$ is prefixed with $P¯$⁠. We can thus calculate $lo(aP¯)=lo(P¯)+∑b<a¯occ(P¯b)=lo(P¯)+∑b<a¯occ(b¯P)$⁠.

It is easy to see that if $(k,k′,s)$ is the SA bi-interval of P, $(k′,k,s)$ is the SA bi-interval of $P¯$⁠, and vice versa. Therefore, we can use the backward extension of $P¯$ to achieve the forward extension of P. Algorithm 3 gives the details. It simplifies the original formulation (Li 2012).

2.6 Finding supermaximal exact matches

An exact match between strings T and P is a 3-tuple (t, p, l) such that $T[t,t+l)=P[p,p+l)$⁠. A maximal exact match (MEM) is an exact match that cannot be extended in either direction. A MEM may be contained in another MEM on the pattern string P. For example, suppose $T=GACCTCCG$ and $P=ACCT$⁠. MEM (5, 1, 2) is contained in MEM (1, 0, 4) on the pattern string. A supermaximal exact match (SMEM) is a MEM that is not contained in other MEMs on the pattern string. In the example above, only (1, 0, 4) is a SMEM. There are usually much fewer SMEMs than MEMs. Gagie (2024) recently proposed a new algorithm to find long SMEMs (Algorithm 4) that is faster than our earlier algorithm (Li 2012) especially when there are many short SMEMs to skip. Both algorithms can also find SMEMs occurring at least $smin$ times (Tatarnikov et al. 2023).

2.7 Finding inexact matches with BWA-SW

The prefix trie of T is a tree that encodes all the prefixes of T (Fig. 1b). Each edge in the tree is labeled with a symbol in Σ. A path from a node to the root spells a substring of T. We can label a node with the SA interval of the string from the node to the root. When we know the label of a node, we can find the label of its child using backward search. We can thus simulate the top-down traversal of the prefix trie (Lam et al. 2008).

As is shown in Fig. 1b, different nodes in the prefix trie may have the same label. If we merge nodes with the same label (Fig. 1c), we will get a prefix DAWG (Blumer et al. 1983). For $|T|>1$⁠, the DAWG has at most $2|T|−1$ nodes (Blumer et al. 1984). Each node uniquely corresponds to an important SA interval of T.

On real data, G_T may be too large to store explicitly. Ropebwt3 instead explicitly stores G_P only and traverses it in the topological order (Algorithm 5). At a node $u∈V(GP)$⁠, we use a hash table to keep ${v∈V(GT):Huv>0}$⁠. This algorithm is exact in that it guarantees to find the best alignment. In practice, however, a large number of $v∈V(GT)$ may be aligned to u with H_uv > 0. It is slow and memory demanding to keep track of all cells (u, v) with positive scores when P is long. Similar to our earlier work, we only store top W cells (25 by default) at each u. This heuristic is akin to dynamic banding for linear sequences (Suzuki and Kasahara 2018).

2.8 Estimating local haplotype diversity

BWA-SW with the banding heuristic may miss the best matching haplotype especially given an index consisting of similar haplotypes. When the suffix on the best full-length alignment has a lot more mismatches than the suffix on suboptimal alignments, the best alignment may have moved out of the band early in the iteration and thus get missed. Nevertheless, a long read sequenced from a new sample may be the recombinant of two genomes in the index. We often do not seek the best alignment of the long read to a single haplotype. We are instead more interested in the collection of haplotypes a query sequence can be aligned to even if they do not lead to the best alignment. This now becomes possible as BWA-SW explores suboptimal alignments to multiple haplotypes.

More exactly, we perform semi-global sequence-to-DAWG alignment (ie requiring the query sequence to be aligned from end to end) by applying BWA-SW to the graph representing the linear query sequence P, from the end to the start. We can find the set of matching haplotypes $M≜{v∈V(GT):Hu0v>0}$ where u₀ represents the start of P. $M$ may include suboptimal haplotypes caused by small variants as well as the optimal one.

When counting hits, we want to ignore the second suboptimal alignment. In theory, we can identify this case by comparing the position of each alignment. This procedure is slow as the locate operation is costly. We instead leverage the bi-directionality to identify redundancy heuristically. Suppose P can be aligned to SA bi-intervals $(k1,k′1,s1)$ and $(k2,k′2,s2)$ and the alignment score to the first interval is higher. We filter out the second interval if $[k2,k2+s2)⊂[k1,k1+s1)$ or $[k′2,k′2+s2)⊂[k′1,k′1+s1)$⁠. This strategy does not avoid all overcounting but it works well on multiple real examples we have closely inspected.

In practice, we may apply this algorithm to the flanking sequence of a variant or to sliding k-mers of a long query sequence to enumerate possible local haplotypes and estimate their frequencies in the index. It is a new query type that is biologically meaningful.

3 Results

3.1 Performance on index construction

We evaluated the performance of BWT construction on 100 haplotype-resolved human assemblies collected in Li et al. (2024). As we included both strands (Section 2.5), each BWT construction algorithm took about 600 billion bases as input (Table 2). grlBWT (commit 5b6d26a; Díaz-Domínguez and Navarro 2023) is the fastest algorithm (Table 3) at the cost of ∼2 terabytes of working disk space including decompressed sequences. Ropebwt3 took 21 h from input sequences in gzip’d FASTA to the final index, of which 7.7 h was spent on libsais. It does not use working disk space and can append new sequences to an existing BWT. However, hardcoded for the DNA alphabet, ropebwt3 does not work with more general alphabets. pfp-thresholds (Rossi et al. 2022) used more memory than the input sequences. It may be impractical with increased sample size.

Both MONI (v0.2.1; Rossi et al. 2022) and Movi (v1.1.0; Zakeri et al. 2024) use pfp-thresholds for building BWT. They generate additional data structures on top of BWT. Time spent on these additional steps were not counted. The tested version of Movi used more than one terabyte of memory to construct the final index. The Movi developer kindly provided the Movi index for the evaluation of query performance in the next section.

We also indexed the same dataset with k-mer based tools. In the lossless mode, metagraph (v0.3.6; Karasikov et al. 2024) indexed the 100 human genomes in 17 h (Table 3). Fulgor (v3.0.0; Fan et al. 2024) is by far the fastest. However, lossy in nature, it is not directly comparable to the rest of the tools in the table.

To test scalability, we indexed a larger dataset consisting of 320 human genomes recently released by the Human Pangenome Reference Consortium (Liao et al. 2023). Because these assemblies are more contiguous than human100 (Table 2), m₂ in Algorithm 1 is smaller, which reduces the multi-threading efficiency of ropebwt3. We can alleviate this by executing BWT merge and partial BWT construction at the same time at the cost of higher memory footprint. On human320, grlBWT is 2–4 times as fast but uses more memory (Table 3). On the largest CommonBacteria dataset (Table 2), ropebwt3 took less than a month to construct the double-strand BWT with 14.66 trillion symbols. The final index in the fermi format is 27.6 GB in size. The “dustbin” and “unknown” categories in AllTheBacteria consist of many unique sequences which are not compressed well. Indexing the entire AllTheBacteria with ropebwt3 will probably take 100–200 GB memory at the peak.

3.2 Query performance

We queried 100–200 Mb human short reads (SR+), human long reads (LR+) and non-human long reads (LR–) against the human pangenome indices constructed above (Table 4). It is important to note that no two tools support exactly the same type of query, but the comparison can still give a hint about the relative performance.

Among the three BWT-based tools, ropebwt3 is slower than Movi but faster than MONI on finding exact matches. Movi finds pseudo-matching length (PML) which is not intended to be the longest exact match. PML corresponds to the longest MEM for only 195 out of one million short reads, and none for long reads. The longest PML of each read is on average 27% shorter than the longest exact match for LR+ and 22% shorter for SR+. Nonetheless, we believe it is possible to implement SMEM finding based on the Movi data structure with minor performance overhead.

MONI and ropebwt3 can also find inexact matches. Implementing r-index, MONI can relatively cheaply locate one SMEM. It leverages this property to extract the genomic sequence around the longest SMEM and performs Smith-Waterman extension. MONI extension is faster than BWA-SW on SR+ because it does not need to inspect suboptimal hits. This feature apparently focuses on short reads only. On LR+, MONI extension fails to extend to the ends of reads and generate incorrect output for the majority of reads.

As to k-mer indices, fulgor outputs the labels of genomes that each read have enough 31-mer matches to. In its current form, such output is not useful for pangenome analysis as we know most of reads can be mapped to all genomes. Metagraph can additionally output the contig name of each match. However, when most k-mers are present in the index, metagraph is impractically slow.

3.3 Identifying novel sequences

As a biological application, we used the pangenome index to identify novel sequences in reads that are absent from other human genomes. For this, we downloaded the PacBio HiFi reads for tumor sample COLO829 (https://downloads.pacbcloud.com/public/revio/2023Q2/COLO829/COLO829/), mapped them to human100 (Table 2), which includes T2T-CHM13 (Nurk et al. 2022), and extracted $≥$1 kb regions on reads that are not covered by SMEMs of 51 bp or longer. We found 95 kb sequences in 43 reads. These sequences could not be assembled. NCBI BLAST suggested multiple weak hits to cow genomes. We could not identify the source of these sequences but there were few of them anyway.

When we mapped the COLO829 reads to T2T-CHM13 only and applied the same procedure, we found 55.9 Mb of “novel” sequences in 25 365 reads. The much larger number is caused by regions with dense SNPs that prevented long exact matches. Counterintuitively, mapping these reads to T2T-CHM13 with minimap2 (Li 2018) resulted in more “novel” sequences at 78.6 Mb, probably because minimap2 ignores seeds occurring thousands of times in T2T-CHM13 and may miss more alignments. Capable of finding SMEMs at the pangenome scale, ropebwt3 is more effective for identifying known sequences. It is also 16% faster and uses less memory than full minimap2 alignment against a single genome.

3.4 Haplotype diversity around C4 genes

The reference human genome GRCh38 has two paralogous C4 genes, C4A and C4B, and they may have copy-number changes (Sekar et al. 2016). In both cases, the exon 26 harbors the functional domain. We extracted the exon 26 from both genes from GRCh38. They differ at six mismatches over 157 bp. We aligned both 157 bp segments, one from each gene, to the human pangenome. 105 haplotypes have the same C4A exon 26 sequence, one haplotype has a mismatch at offset 128 and four haplotypes have a mismatch at 54. In case of C4B, 60 haplotypes have the same reference C4B sequence and 23 also have a mismatch at offset 54. This means among the 100 human haplotypes, there are 110 C4A gene copies and 83 C4B copies. If we increase the bandwidth from the default 25 to 200, BWA-SW will be able to align C4A exon 26 to C4B genes and output all five hits for each sequence.

Figure 3 shows the local haplotype diversity across the entire C4A gene spanning ∼20.6 kb on GRCh38. We can see most regions have 193 copies, except a ∼6.4 kb HERV insertion that separates long and short forms (Sekar et al. 2016). The dip around exon 26 is caused by the C4A–C4B difference. We could only see these alternative haplotypes with BWA-SW which reports suboptimal hits.

4 Discussions

Ropebwt3 is a fast tool for BWT construction and sequence search for redundant DNA sequences. Generating BWT purely in memory and supporting incremental build, ropebwt3 is convenient and practical for BWT construction at large scale. It is the only algorithm that can construct the BWT of 320 human genomes from a 2.5 GB AGC archive (Deorowicz et al. 2023) without staging all decompressed data in memory or on disk. It provides the fastest algorithm so far for finding supermaximal exact matches and can report inexact hits as well. Ropebwt3 demonstrates that BWT-based data structures are scalable to terabases of pangenome data.

Ropebwt3 implements an FM-index to locate SMEMs or local hits. Although the standard r-index is faster than an FM-index of the same size, it imposes a fixed sampling rate: two suffix array values per run. The BWT of CommonBacteria has 14.6 trillion bases and 17.6 billion runs. An r-index is likely to take >200 GB, while an FM-index sampled at one position per 8192 bp takes 17.5GB in ropebwt3. Subsampled r-index (sr-index; Cobas et al. 2024) is probably the better solution, which we will explore in future.

This article has not evaluated several recent BWT-based tools including r-index-f (Brown et al. 2022), block_RLBWT (Díaz-Domínguez et al. 2023), Move-r (Bertram et al. 2024) and b-move (Depuydt et al. 2024). Although they have not been tested on large datasets like human pangenome and they do not report SMEMs, some of their data structures are probably more efficient than the ones we use. We may integrate these methods into ropebwt3 in future as well.

We often use pangenome graphs to analyze multiple similar genomes (Liao et al. 2023). These graphs are built from the multiple sequence alignment through complex procedures involving many parameters (Li et al. 2020, Garrison et al. 2024, Hickey et al. 2024). It is challenging to understand if the graph topology is biologically meaningful especially given that we do not know the correct alignment between two genomes, let alone multiple ones. Complement to graph-based data structures, BWT-based algorithms are often exact with no heuristics or parameters but they tend to support limited query types. For example, we cannot project the alignment to a designated reference genome. What additional query types we can achieve will be of great interest to the comprehensive pangenome analysis in future.

Acknowledgements

We are grateful to Travis Gagie for pointing us to his long MEM finding algorithm, to Ilya Grebnov for adding the support of 16-bit alphabet which helps to accelerate ropebwt3, to Mohsen Zakeri for providing the Movi index, to Massimiliano Rossi for explaining the MONI algorithm, and to Giulio Pibiri and Rob Patro for trouble-shooting compilation issues with fulgor.

Conflict of interest

None declared.

Funding

This work was supported by National Institute of Health grant [R01HG010040 and U01HG010961 to H.L.].

Data availability

The ropebwt3 source code is available at https://github.com/lh3/ropebwt3. Prebuilt ropebwt3 indices can be obtained from https://doi.org/10.5281/zenodo.11533210 and https://doi.org/10.5281/zenodo.13955431.

References