b-move: faster lossless approximate pattern matching in a run-length compressed index

Abstract

Background: Due to the increasing availability of high-quality genome sequences, pan-genomes are gradually replacing single consensus reference genomes in many bioinformatics pipelines to better capture genetic diversity. Traditional bioinformatics tools using the FM-index face memory limitations with such large genome collections. Recent advancements in run-length compressed indices like Gagie et al.'s r-index and Nishimoto and Tabei's move structure, alleviate memory constraints but focus primarily on backward search for MEM-finding. Arakawa et al.'s br-index initiates complete approximate pattern matching using bidirectional search in run-length compressed space, but with significant computational overhead due to complex memory access patterns.

Results: We introduce b-move, a novel bidirectional extension of the move structure, enabling fast, cache-efficient, lossless approximate pattern matching in run-length compressed space. It achieves bidirectional character extensions up to 7 times faster than the br-index, closing the performance gap with FM-index-based alternatives. For locating occurrences, b-move performs $\varphi$ and ${\varphi }^{-1}$ operations up to 7 times faster than the br-index. At the same time, it maintains the favorable memory characteristics of the br-index, for example, all available complete E. coli genomes on NCBI's RefSeq collection can be compiled into a b-move index that fits into the RAM of a typical laptop.

Conclusions: b-move proves practical and scalable for pan-genome indexing and querying. We provide a C++ implementation of b-move, supporting efficient lossless approximate pattern matching including locate functionality, available at https://github.com/biointec/b-move under the AGPL-3.0 license.

Keywords: Approximate pattern matching; Bidirectional search; Cache efficiency; FM-index; Lossless alignment; Move structure; Pan-genomics; r-index.

Algorithm 1

Find $\text{LF}(i)=M_{\text{LF}}.\mathtt{moveStep}(i,j)$ given tuple (i, j), where j is the run index that contains BWT index i.

Algorithm 2

Let [s, e] be the interval over $\text{SA}$ corresponding to P, which we want to extend to cP. $R_s$ and $R_e$ are the runs containing s and e, respectively. Functions walkToNextRun and walkToPreviousRun return the $\text{SA}$ indices $s_c$ and $e_c$ (along with their run indices $R_{s,c}$ and $R_{e,c}$) which indicate the smallest subinterval within [s, e] containing all occurrences of c in the BWT.

Algorithm 3

Extend pattern P, represented by [s, e], corresponding to respectively run index $R_s$ and $R_e$, to cP. The algorithm returns the updated interval $[s^{\prime },e^{\prime }]$, corresponding to respectively run index $R_s^{\prime }$ and $R_e^{\prime }$, for cP.

Algorithm 4

Update intervals $([s,e],R_s,R_e)$ and $([s^{\text{rev}},e^{\text{rev}}],R_s^{\text{rev}},R_e^{\text{rev}})$ corresponding to P, to intervals $([s^{\prime },e^{\prime }],R_s^{\prime },R_e^{\prime })$ and $([{s^{\text{rev}}}^{\prime },{e^{\text{rev}}}^{\prime }],R_s^{\text{rev}},R_e^{\text{rev}})$ that correspond to cP.

Algorithm 5

Update run indices using binary search to find the correct run containing index i. The algorithm uses the move table $M_{\text{LF}}$ upon which it is called.

Fig. 1

Log scale histograms for the number of walking and fast forwarding steps required per successful bidirectional character extension (96,363,328 in total) for finding all SA intervals corresponding to all occurrences up to 3 mismatches. We aligned 100000 Illumina reads of length 151 bp and their reverse complements to the pan-genome composed of 512 human chromosome 19 haplotypes

Fig. 2

Benchmark results for finding all SA intervals corresponding to all occurrences up to a certain number of mismatches (x-axis). The left panel displays the total search time, while the right panel shows the average execution time for the core bidirectional character extension functionality. We aligned 100000 Illumina reads of length 151 bp and their reverse complements to the pan-genome composed of 512 human chromosome 19 haplotypes

Fig. 3

Log scale histogram for the number of fast forwarding steps required per $\varphi$ or ${\varphi }^{-1}$ operation (355,041,072 in total) for finding all occurrences up to 3 mismatches. We aligned 100000 Illumina reads of length 151 bp and their reverse complements to the pan-genome composed of 3264 E. coli strains, using unbalanced $\varphi$ move tables. Using balanced $\varphi$ move tables, the mean number of fast forwarding steps is reduced to 1.49 (not shown)

Fig. 4

Benchmark results for finding all occurrences with up to a certain number of mismatches (x-axis). The left panel presents the average locate time per SA interval, while the right panel shows the average execution time for the core $\varphi$ and ${\varphi }^{-1}$ operations. The experiment aligns 100,000 Illumina reads of length 151 bp and their reverse complements to the pan-genome composed of 3264 E. coli strains

Fig. 5

Benchmark results for approximate pattern matching using the br-index, b-move (with and without $\varphi$ move tables), and Columba (suffix array sparseness of 8). Additionally, we include results for b-move (with $\varphi$ move tables) with reporting functionality, which is not included in the other results. We aligned 100000 Illumina reads of length 151 bp and their reverse complements to pan-genomes for both human chromosome 19 and E. coli, across multiple sizes. We allowed for a maximum error distance of 3 (Hamming for br-index, edit for b-move and Columba)

Fig. 6

Log-scale time-memory comparison for approximate pattern matching using the br-index, b-move (with and without $\varphi$ move tables), and Columba (suffix array sparseness of 8). We show the switch point for each pan-genome, which marks the minimum pan-genome size at which the b-move index becomes smaller than Columba. We aligned 100000 Illumina reads of length 151 bp and their reverse complements to pan-genomes for both human chromosome 19 and E. coli, across multiple sizes. We allowed for a maximum error distance of 3 (Hamming for br-index, edit for b-move and Columba)