Bioinformatics, 39(2), 2023, btad074 https://doi.org/10.1093/bioinformatics/btad074 Advance Access Publication Date: 7 February 2023 Original Paper Sequence analysis Optimal gap-affine alignment in O(s) space Santiago Marco-Sola 1,2,*, Jordan M. Eizenga 3, Andrea Guarracino4,5, Benedict Paten 3, Erik Garrison 5 and Miquel Moreto1,6 1 Computer Sciences Department, Barcelona Supercomputing Center, Barcelona 08034, Spain, 2Departament d’Arquitectura de Computadors i Sistemes Operatius, Universitat Autònoma de Barcelona, Barcelona 08193, Spain, 3Genomics Institute, University of California Santa Cruz, Santa Cruz, CA 95064, USA, 4Genomics Research Centre, Human Technopole, Milan 20157, Italy, 5Department of Genetics, Genomics and Informatics, University of Tennessee Health Science Center, Memphis, TN 38163, USA and 6Departament d’Arquitectura de Computadors, Universitat Politècnica de Catalunya, Barcelona 08034, Spain *To whom correspondence should be addressed. Associate Editor: Pier Luigi Martelli Received on August 17, 2022; revised on January 2, 2023; editorial decision on January 30, 2023 Abstract Motivation: Pairwise sequence alignment remains a fundamental problem in computational biology and bioinfor- matics. Recent advances in genomics and sequencing technologies demand faster and scalable algorithms that can cope with the ever-increasing sequence lengths. Classical pairwise alignment algorithms based on dynamic pro- gramming are strongly limited by quadratic requirements in time and memory. The recently proposed wavefront alignment algorithm (WFA) introduced an efficient algorithm to perform exact gap-affine alignment in OðnsÞ time, where s is the optimal score and n is the sequence length. Notwithstanding these bounds, WFA’s Oðs 2 Þ memory requirements become computationally impractical for genome-scale alignments, leading to a need for further improvement. Results: In this article, we present the bidirectional WFA algorithm, the first gap-affine algorithm capable of comput- ing optimal alignments in OðsÞ memory while retaining WFA’s time complexity of OðnsÞ. As a result, this work improves the lowest known memory bound OðnÞ to compute gap-affine alignments. In practice, our implementation never requires more than a few hundred MBs aligning noisy Oxford Nanopore Technologies reads up to 1 Mbp long while maintaining competitive execution times. Availability and implementation: All code is publicly available at https://github.com/smarco/BiWFA-paper. Contact: santiagomsola@gmail.com Supplementary information: Supplementary data are available at Bioinformatics online. 1 Introduction these methods compute the optimal alignment score for progressive- Pairwise sequence alignment provides a parsimonious transformation ly longer prefixes of q and t, which correspond to the cells of the DP of one string into another. From this transformation, we can under- matrix. Thus, an optimal alignment can then be read out by tracing stand the relationship between pairs of sequences. Because similarities the recurrence back through the matrix. and differences between biosequences (DNA, RNA, protein) relate to Selecting a suitable alignment score function is essential to obtain variation in function and evolutionary history of living things, pair- biologically meaningful alignments, as it determines the characteristics wise sequence alignment algorithms are a core part of many essential of optimal alignments. In effect, the alignment score function encodes bioinformatics methods in read mapping (Li, 2013; Marco-Sola et al., prior expectations about the probability of certain kinds of sequence 2012), genome assembly (Koren et al., 2017; Simpson et al., 2009), differences. It has been observed that, in many contexts, insertions variant calling (Garrison and Marth, 2012; McKenna et al., 2010; and deletions are non-uniformly distributed; they are infrequent but Rodrıguez-Martın et al., 2017) and many others (Durbin et al., 1998; tend to be adjacent so that they form extended gaps with a long-tailed Jones et al., 2004). Its importance has motivated the research and de- length distribution. This motivated the development of gap-affine velopment of multiple solutions over the past 50 years. models in which the penalty of starting a new gap is larger than that Classical approaches to derive alignments involve the application of extending a gap (Gotoh, 1982). Crucially, gap-affine penalties can of dynamic programming (DP) techniques. These methods require be implemented efficiently using additional DP matrices. computing a matrix whose dimensions correspond to the lengths of Problematically, the efficiency of classical gap-affine DP-based the query q and target t sequences. Using DP recurrence relations, methods is constrained by their quadratic requirements in time and C The Author(s) 2023. Published by Oxford University Press. V 1 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unre- stricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 2 S.Marco-Sola et al. memory with respect to the lengths of the sequence pair. character. We will use ðx; o; eÞ to denote the gap-affine penalties. A Consequently, multiple optimizations have been proposed over the mismatch costs x, and a gap of length l costs o þ l  e. We assume years. Notable examples include bit-parallel techniques (Loving that x > 0 and e > 0, and further that all of the score parameters are et al., 2014), data-layout transformations to exploit SIMD instruc- constants. tions (Farrar, 2007; Rognes and Seeberg, 2000; Wozniak, 1997), Basically, WFA computes partial optimal alignments of increas- difference encoding of the DP matrix (Suzuki and Kasahara, 2018), ing score until an alignment with score s reaches coordinate ðn; mÞ among other methods (Altschul et al., 1990; Kiełbasa et al., 2011; of the DP matrix. In this way, the algorithm determines that s is the Xia et al., 2021; Zhao et al., 2013). Nonetheless, all these exact minimal alignment score. Moreover, it can derive the optimal align- methods retain the quadratic requirements of the original DP algo- ment by tracing back the partial alignments that led to score s at rithm and therefore struggle to scale when aligning long sequences. ðn; mÞ. In many cases, when two sequences are homologous, the majority Let Ms;k , X s;k , I s;k and Ds;k denote the offset within diagonal k of possible alignments are largely sub-optimal, having a substantially in the DP-matrix to the farthest-reaching (f.r.) cell that has score s worse score than the optimal one. For this reason, heuristic methods and ends with a match, mismatch, insertion or deletion, respectively. are usually employed to find candidate alignment regions when the In general, we denote by wavefront the tuple of offsets for a given cost of exact algorithms becomes impractical. Most notable score W s ¼ ðMs ; X s ; I s ; Ds Þ. We refer to the four elements in this approaches use adaptive band methods (Suzuki and Kasahara, 2017) tuple as its components, and we associate a corresponding sentinel or pruning strategies [e.g. X-drop (Zhang et al., 2000) and Z-drop value to specify each component: c 2 fM; X; I; Dg. (Li, 2018)] to avoid the computation of alignments extremely unlikely In Marco-Sola et al. (2021), we proved that the f.r. points of W s to be optimal. These heuristic methods have been implemented within can be computed using previous wavefronts W soe , W se and many widely used tools (Altschul et al., 1990; Li, 2018). W sx , using Equation 1 where LCPðv; wÞ is the length of longest Recently, we proposed the wavefront alignment algorithm common prefix between substrings v and w. The base case for this (WFA) (Marco-Sola et al., 2021) to compute the exact alignment be- recursion is given by X 0;0 ¼ 0. tween two sequences using gap-affine penalties. WFA reformulates the alignment problem to compute the longest-possible alignments I s;k ¼ maxfMsoe;k1 þ 1; I se;k1 þ 1g of increasing score until the optimal alignment is found. Notably, Ds;k ¼ maxfMsoe;kþ1 ; Dse;kþ1 g (1) WFA takes advantage of homologous regions between sequences to X s;k ¼ maxfMsx;k þ 1; I s;k ; Ds;k g accelerate alignment’s computation. As a result, WFA computes op- Ms;k ¼ X s;k þ LCPðq½X s;k  k; n  1; t½X s;k ; m  1Þ; timal gap-affine alignments in OðnsÞ time and Oðs2 Þ memory, where Equation 1 shows that the computation of a given wavefront n is the sequence length and s the optimal alignment score. Being an depends only on the previous p ¼ maxfx; o þ eg wavefronts. We exact algorithm, WFA provides the same guarantee for optimality as refer to p as the wavefront scope or, in other words, the maximum classical algorithms (Gotoh, 1982; Needleman and Wunsch, 1970; score increase between partial alignments. Moreover, note that X s;k Smith and Waterman, 1981), but does away with the quadratic does not need to be explicitly stored as its values can be inferred requirements in time. using Ms;k , I s;k and Ds;k . WFA unlocked the path for optimal alignment methods capable In the worst case, WFAPrequires computing s wavefronts of of scaling to long sequences. Nevertheless, the Oðs2 Þ memory increasing length, totalling s 2 i¼0 ð1 þ 2iÞ ¼ Oðs Þ cells. Moreover, requirements quickly become the limiting factor when aligning suffi- the LCP must be computed once for each cell. However, within a di- ciently long or noisy sequences (Eizenga and Paten, 2022). As it hap- agonal, the total number of offset increments cannot exceed the pens, WFA’s memory requirements can be impractical when length of the sequences. Hence, WFA requires OðnsÞ time and Oðs2 Þ aligning through large structural variations or highly divergent gen- memory in the worst case (Marco-Sola et al., 2021). Since s  pn, ome regions. Given that we use alignment to understand variation, the OðnsÞ factor of the execution time, due to the LCP, dominates these are some contexts in which optimal alignment could be most over the Oðs2 Þ factor in the worst case. However, in practice, the useful, but its memory requirements make it prohibitive. time is often closer to Oðs2 þ nÞ. This is because spurious matches To address this problem, this article presents the first gap-affine between high-entropy sequences are short in expectation. alignment algorithm to compute the optimal alignment in OðnsÞ time Accordingly, the LCP computations often finish after performing and OðsÞ memory (excluding the storage of the input sequences). Our only a few character comparisons, except along the optimal align- method, the bidirectional WFA algorithm (BiWFA), computes the ment in which OðnÞ comparisons are required. WFA alignment of two sequences in the forward and reverse direction until they meet. Using two wavefronts of OðsÞ memory, we demon- strate how to find the optimal breakpoint of the alignment at score 2.2 Bidirectional wavefront alignment algorithm s=2 and proceed recursively to solve the complete alignment in The core idea of the BiWFA algorithm is to perform WFA simultan- OðnsÞ time. To our knowledge, this work improves the lowest known eously in both directions on the strings: from start to end (i.e. for- memory bound to compute gap-affine alignments OðnÞ (Myers and ward) and from end to start (i.e. reverse). Each direction will only Miller, 1988) to OðsÞ, while retaining the time complexity of the ori- retain p wavefronts in memory. This is insufficient to perform a full ginal WFA OðnsÞ. Furthermore, our experimental results demonstrate traceback. However, when they ‘meet’ in the middle, we can infer a that the BiWFA delivers comparable, or even better, performance breakpoint in the alignment that divides the optimal score roughly than the original WFA, outperforming other state-of-the-art tools in half. Then, we can apply the same procedure on the two sides of while using a minimal amount of memory. the breakpoint recursively. We will show that this results in only a The rest of the article is structured as follows. Section 2 presents constant-factor slowdown. This technique was previously employed the definitions, algorithms and formal proofs supporting BiWFA. to a similar end with the Myers OðNDÞ difference algorithm Section 3 shows the experimental evaluation of our method, com- (Myers, 1986). paring it against other state-of-the-art tools and libraries. Lastly, Figure 1 presents a graphical example of BiWFA computing a Section 4 presents a discussion on the BiWFA method and summa- breakpoint in the optimal alignment between two sequences. The rizes the contributions and impact of this work. figure shows the DP cells computed by the forward and reverse wavefronts. Alignments in both directions progress until they over- lap on cell ð4; 4Þ with score 8 þ 8 ¼ 16 corresponding to the optimal 2 Materials and methods alignment (sopt ¼ 16). First, let us define the WFA equations for the forward and re- 2.1 Wavefront alignment algorithm verse alignment directions. The recursions for the forward direction Let the query q ¼ q0 q1 . . . qn1 and the text t ¼ t0 t1 . . . tm1 be are equivalent to those of the standard WFA presented above strings of length n and m, respectively. Likewise, let v½i; j ¼ (Equation 1). However, to highlight the distinction, we will denote ! ! ! ! ! vi viþ1 . . . vj denote a substring of any string v from the ith to the jth them W s ¼ ð I s ; D s ; X s ; M s Þ. The recursions for the reverse Optimal gap-affine alignment in O(s) space 3  ! previously computed reverse wavefronts W sr ...sr p . First, if W sf belongs to a score-balanced breakpoint (with jsf  sr j  p), it is suf-  ficient to check for overlaps against W sr and the previous p  1 re- ! verse wavefronts. Second, for every diagonal k in wavefront W sf , Algorithm 2 checks of overlaps in all wavefront components. This way, the algorithm keeps track of the overlap with the minimum score detected so far. Last, note that overlaps on I and D compo- nents account twice for the gap-open score o. Hence, the score from overlaps at indel components has to be decreased by o. Fig. 1. Example of BiWFA aligning q ¼ “TCTAGCG” against t ¼ “TGGAAAG” under the penalties ðx ¼ 4; o ¼ 5; e ¼ 1Þ direction are very similar (Equation 2), using X 0;mn ¼ m as the base case and LCSðv; wÞ to denote the length of the longest common suffix of v and w. Note that the same argument used in Marco-Sola et al. (2021) applies to the reverse recursions to prove that they are f.r. in the reverse direction.  I s;k ¼ minf M soe;kþ1  1; I se;kþ1  1g  D s;k ¼ minf M soe;k1 ; D se;k1 g  (2) X s;k ¼ minf M sx;k  1; I s;k ; D s;k g  M s;k ¼ X s;k  LCSðq½0; X s;k  k  1; t½0; X s;k  1Þ Algorithm 1 presents the BiWFA algorithm to compute a break- point in the optimal alignment at s=2. Using forward and reverse wavefronts, the algorithm proceeds by alternatingly computing !  !  forward and reverse alignments (i.e. W 1 ; W 1 ; W 2 ; W 2 ; . . .). To this end, BiWFA relies on the operators WF_NEXT() and WF_EXTEND() from the standard WFA [see Marco-Sola et al. (2021)] to compute successive wavefronts using Equaions 1 and 2. The process is halted after their offsets overlap to compute the pos- ition of a breakpoint in the optimal alignment. This algorithm iter- ates until it is guaranteed that the optimal breakpoint has been found. However, there are some technical details involving the de- tection of overlaps and the computation of the optimal breakpoint, which we cover in Sections 2.3 and 2.4. 2.3 Finding a score-balanced breakpoint in the optimal alignment The first technical detail involved in finding an alignment break- point between the two directions is that it is often not possible to split an alignment into an equally scoring prefix and suffix. In gen- eral, two prefixes of the optimal alignment that differ by one charac- ter can have scores that differ by as much as p. Accordingly, we will demand a weaker notion of balance. If sf and sr are the forward and reverse scores, respectively, we will aim to have jsf  sr j  p. The second technical detail is that the optimal score is not always the sum of the two scores. This occurs because the forward iteration In practice, Algorithm 1 can avoid most calls to incurs the gap open penalty o at the beginning of gaps, but the re- WF_OVERLAP(). An efficient implementation can keep track of verse incurs it at the end of gaps (or rather, at the beginning in the the farthest reached antidiagonal by each wavefront. If the most reverse direction). Thus, if the two directions meet in a gap, then we advanced antidiagonal reached by the forward and reverse wave- have sopt ¼ sf þ sr  o rather than sopt ¼ sf þ sr , where sopt is the op- fronts do not overlap (ANTIDIAG() on Algorithm 1), it follows timal alignment score. that no offsets from any diagonal can overlap, rendering the call to The final technical detail is that offsets of the two directions may WF_OVERLAP() unnecessary. not precisely meet. WFA proceeds by greedily taking matches in both directions. This makes it possible for the two directions to shoot past each other without actually meeting. It turns out that it is 2.4 Correctness of the breakpoint detection sufficient to detect that such an overshoot has occurred, as will be The correctness of the Algorithm 1 stems from the following lemma. shown in Section 2.4. In Algorithm 2, we reconcile these three difficulties. Without loss Lemma 2.1. The optimal alignment score sopt  s if and only if ! of generality, we assume that a forward wavefront W sf has been there exist sf , sr and k such that jsf  sr j  p and at least one of the computed (Algorithm 1), and we want to detect overlaps against the following is true: 4 S.Marco-Sola et al. !  1. sf þ sr ¼ s and M k;sf  M k;sr 2.6 BiWFA uses OðsÞ space and Oððm þ nÞsÞ time ! The memory complexity of Algorithm 1 is relatively simple to char- 2. sf þ sr ¼ s þ o and I k;sf  I k;sr acterize. The range of diagonal values k increases by at most two ! every time s is incremented, and each forward and reverse search 3. sf þ sr ¼ s þ o and D k;sf  D k;sr  only needs to store the last p wavefronts. Thus, the memory use is and further, M k;sr (resp. I k;sr , D k;sr ) is included in the traceback of proportional to the optimal alignment score, OðsÞ, excluding the an alignment with score at most s if the first (resp. second, third) storage of the input sequences. Also, note that the output alignment condition is true. only requires storing the position ði; jÞ for the mismatches, insertions Proof. and deletions (matches can be inferred from the gaps). Concerning Algorithm 3, data structures are discarded before entering a recur- See Supplementary Material. sive call. Therefore, the maximum memory use occurs in the outer- most call, in which s is the optimal score of the full alignment. The time complexity is more complicated to analyze. Our proof follows similar arguments as those from Myers (1986). This lemma implies that the minimum value s for which the ‘only if’ condition holds is the optimal score. Moreover, if the first of Theorem 2.2. BiWFA’s time complexity is Oððm þ nÞsÞ, being n and the three conditions is found to hold for some values sf and sr , m the sequences’ length and s the optimal alignment score. then sopt  sf þ sr þ o. Therefore, Algorithm 1 is guaranteed to find part of a minimum-scoring alignment based on the following Proof. features: Let ‘ ¼ m þ n, and let Tð‘; sÞ be BiWFA’s execution time with score s. A call to BiWFA can result in two recursive calls. Let ‘f and ‘r be • Algorithm 2 checks a window of p score values on each iteration. the combined length of the sequences in the two calls, and similarly let sf and sr be the two alignment scores. Following Lemma 2.1, we • Algorithm 1 iterates through alternatingly increasing values of sf know that these variables obey the following inequalities: and sr , detecting breakpoints with scores of at least sf þ sr  o  p þ 1 in each iteration. • ‘f þ ‘r  ‘ • After finding some sf and sr that satisfy the overlap condition, Algorithm 1 continues for additional iterations until it is no longer • sf þ sr  s possible to find a lower score. • jsf  sr j  p Because each direction of WFA executes in Oðs‘Þ time (Marco-Sola et al., 2021), we can choose a constant c1 large enough that the fol- 2.5 Combining breakpoints into an alignment lowing inequality holds for all s > 3p: Algorithm 3 shows how to use BiWFA to recursively split align- ments into smaller subproblems until the remaining alignment can be trivially solved. Tð‘; sÞ  c1 s‘ þ Tð‘f ; sf Þ þ Tð‘r ; sr Þ (3) Note that a breakpoint computed by BiWFA can be found on the I or D components. Thus, those alignments that connect with this We can also choose a constant c2 large enough that for all s  3p breakpoint have to start or end at the given component. This way, Algorithm 3 considers the starting and ending component of each Tð‘; sÞ  c2 ‘ (4) alignment, and forces the underlying WFAs to use different initial conditions depending on the alignment starting at the M (X 0;0 ¼ 0), I (I 0;0 ¼ 0) or D component (D0;0 ¼ 0). A similar argument applies This follows because the recursion depth depends only on s, which to the ending conditions of each alignment ending at the M we have given an upper bound. Therefore, this term includes a (Ms;mn ¼ m), I (I s;mn ¼ m) or D component (Ds;mn ¼ m). bounded number of calls that all have linear dependence on ‘. Optimal gap-affine alignment in O(s) space 5 Next, we argue that Tð‘; sÞ  3c1 s‘ þ c2 ‘ by induction on s. The the executions. Therefore these methods were discarded from the base cases for s ¼ 0; 1; . . . ; 3p follow trivially from the latter of the evaluation. previous inequalities. Assume then that s > 3p and the induction hy- All the presented methods have been configured to generate glo- pothesis holds for 0; 1; . . . ; s  1. Note that we then have bal alignments. These algorithms are grouped in two categories: sf ; sr  2s=3, else either jsf  sr j > p or sf þ sr > s. Thus, ‘Gap-affine Exact’ for exact algorithms that use gap-affine penalties (i.e. BiWFA, WFA and its low-memory modes, wfalm and its low- Tð‘; sÞ  c1 s‘ þ ð3c1 ð2s=3Þ‘f þ c2 ‘f Þ þ ð3c1 ð2s=3Þ‘r þ c2 ‘r Þ memory modes and KSW2-Z2), and ‘Others’ for methods that use (5)  3c1 s‘ þ c2 ‘ simpler penalty models or can only compute the alignment score This proves the claim. (i.e. Edlib and BitPal). For the evaluation, we use simulated and real datasets. For the simulated datasets, we simulate several datasets of various sequence lengths (i.e. 100K, 500K, 1M and 2M bases) and different error rate 3 Results (i.e. e ¼ 10% and 20%) randomly generated. Regarding the evalu- We implement the BiWFA algorithm described in this work in C. ation with real datasets, we use a first set of sequences generated by The code and the scripts required to reproduce the experimental the Human Pangenome Reference Consortium (Miga and Wang, results presented in this section are publicly available and can be 2021), consisting of long reads sequenced using Oxford Nanopore found at https://github.com/smarco/BiWFA-paper. Moreover, the Technologies (ONT), PromethION platform, with an average error code has been integrated into the WFA2-lib alignment library (as rate of 5–10%. The sequences are derived from the human cell line ultralow memory mode) at https://github.com/smarco/WFA2-lib. HG002, subset to chromosome 12 and restricted to those at least 10 kbp long, for a total number of 1312 sequence pairs of average 3.1 Experimental setup length equal to 172 kbp (maximum 306 kbp). In addition, we use a We evaluate the performance of our BiWFA implementation com- second dataset comprising ONT MinION reads from Bowden et al. pared to the state-of-the-art and other high-performance sequence (2019), with an average error rate of 5% and restricted to those at alignment libraries. We select the original WFA (Marco-Sola et al., least 500 kbp long, for a total number of 48 sequence pairs of aver- 2021) (wfa-high) and its new low-memory modes (wfa-med and age length equal to 630 kbp (maximum 1 Mbp). All the executions are performed single-thread on a node running wfa-low) implemented in WFA2-lib (https://github.com/smarco/ CentOS Linux (release 8.1.1911) equipped with an AMD EPYC WFA2-lib). Also, we select the efficient wfalm (Eizenga and Paten, 7742 CPU and 1 TB of DRAM (distributed in 16 dimms  64 GiB 2022) (wfalm) and its low-memory modes (wfalm-low and wfalm- @3200 MHz). rec). Moreover, we include the highly optimized KSW2-Z2 (ksw2_extz2_sse), from the KSW2 library (Li, 2018; Suzuki and Kasahara, 2018), as the best representative of DP-based methods 3.2 Evaluation on simulated data due to its exceptional performance and widely usage within bio- Table 1 shows the performance results (i.e. execution time and informatics tools. In addition, we include the Edlib (Sosic and  Sikic, memory) for the different methods using simulated datasets. 2017) and BitPal (Loving et al., 2014) libraries, which implement Overall, the results show that BiWFA is faster and uses less memory bit-parallel alignment strategies for edit and non-unitary penalties than all other methods in the ‘Gap-affine Exact’ category. In par- (i.e. gap-linear), respectively. Although they solve a considerably ticular, BiWFA requires 32  1000 less memory than KSW2-Z2, easier problem (i.e. Edlib is restricted to edit-alignments and BitPal while being 1:4  4:7 faster. Compared to original WFA-based only computes the alignment score), and thus are not directly com- methods (i.e. WFA-high and wfalm), BiWFA uses 9  9620 less parable, we included them in the evaluation to provide a perform- memory, being up to 4:4 faster. Similarly, BiWFA outperforms the ance upper bound. other memory-efficient WFA-based methods (i.e. WFA-med, WFA- We considered including other popular methods like those low, wfalm-low and wfalm-rec), reducing memory requirements implemented in the Parasail (Daily et al., 2015; Daily, 2016; Farrar, down to 438 while being 2:7  26:5 faster. More importantly, 2007; Wozniak, 1997), SeqAn (Rahn et al., 2018) and Gaba (Suzuki most of the pairwise alignment methods evaluated fail to scale and Kasahara, 2018) libraries. However, these libraries were not megabases-long sequences, requiring more memory than available designed to align long and noisy sequences, and failed to complete in the node (i.e. 1TB). As opposed, BiWFA only requires a few Table 1. Time and memory performance of pairwise alignment implementations on simulated data Time (s) Memory (MBs) 10 Kbp 100 Kbp 1 Mbp 2 Mbp 10 Kbp 100 Kbp 1 Mbp 2 Mbp 10% 20% 10% 20% 10% 20% 10% 20% 10% 20% 10% 20% 10% 20% 10% 20% edlib 0.4 0.6 2.5 4.5 17.9 35.3 35.4 69.1 4 4 5 5 13 13 22 23 bitpal 1.3 1.2 12.3 12.3 123.8 123.7 248.0 247.1 4 4 4 6 10 10 15 13 ksw2-extz2 9.8 9.9 96.7 97.5 n/a n/a n/a n/a 193 196 19 081 19 083 n/a n/a n/a n/a WFA-high 2.0 5.7 28.5 84.1 312.9 n/a n/a n/a 128 313 8981 26 667 932 199 n/a n/a n/a WFA-med 6.6 20.2 89.5 272.8 1922.1 3690.1 n/a n/a 35 81 830 1620 42 464 24 874 n/a n/a WFA-low 8.0 24.1 101.4 301.1 4394.9 4857.0 7710.2 9813.1 25 60 554 884 25 321 12 539 52 551 26 067 wfalm 6.4 19.2 90.2 268.7 841.0 n/a n/a n/a 54 148 8968 26 575 898 770 n/a n/a n/a wfalm-low 10.1 30.4 164.1 494.8 1525.1 4418.7 2990.1 8779.4 10 16 443 823 10 435 30 817 36 299 69 312 wfalm-rec 22.3 70.7 447.5 1402.6 5792.7 17 752.9 11 979.6 37 747.7 6 7 43 73 497 904 1064 1787 BiWFA 2.4 6.9 20.8 61.0 218.3 680.4 466.9 1429.0 6 5 19 27 97 180 202 267 BiWFA.score 1.1 3.1 10.2 30.2 112.2 355.8 245.4 750.2 3 4 16 23 97 186 204 256 Note: Execution time (in seconds) and memory (in MBs) required per 1M bases aligned, using different pairwise alignment implementations on simulated data- sets. Executions that failed appear as ‘n/a’. Best performing implementation in the ‘Gap-affine Exact’ category is marked in bold. Although Edlib and BitPal are not directly comparable to the other implementations, we included them in the comparison as a reference. Similarly, we include executions of BiWFA limited to compute the alignment score as ‘BiWFA.score’. 6 S.Marco-Sola et al. hundred MBs of memory. Note that, computing the full alignment even faster, 1:1  1:28 faster on average) computing exact align- (BiWFA) requires a similar amount of memory as computing only ments (not just the score) under the gap-affine model. the alignment score (BiWFA.score). Nonetheless, computing the For completeness, Supplementary Figure S1 shows experimental alignment score is 2 faster than computing the full alignment. results limited to aligning sequences up to 10 Kbps. In this scenario, For completeness, we present a comprehensive experimental BiWFA demonstrates to be one of the fastest implementations, evaluation on a wider range of sequence lengths (i.e. 100 bp, 1 Kbp, requiring less than 10 MB to execute. 10 Kbp, 100 Kbp, 1 Mbp and 2 Mbp) and error rates (0.1%, 1%, 5%, 10%, 20% and 40%) on Supplementary Material (Supplementary Tables S1 and S2). For short sequences (i.e. 4 Discussion 1 Kbp), the results show that BiWFA delivers similar performance As long sequencing technologies improve and high-quality sequence as the original WFA (1:25  2 slower) while reducing the memory assembly decreases in cost, we anticipate that the importance of requirements up to 60. Our experiments indicate that BiWFA pairwise alignment algorithm will continue to increase. To keep up starts outperforming the original WFA when aligning sequences lon- with upcoming improvements in sequencing and genomics, pairwise ger than 30 Kbp. Nevertheless, the exact performance breakpoint alignment algorithms need to face crucial challenges in reducing exe- can vary depending on the error rate, implementation and processor cution time and memory consumption. In this work, we have pre- specifics. Similarly, for smaller error rates (i.e.  1%), WFA-based sented the BiWFA, a gap-affine pairwise alignment algorithm that methods largely outperform other approaches, being 2  3 orders of requires OðnsÞ time and OðsÞ space, being the first algorithm to im- magnitude faster than other methods like KSW2-Z2, bitpal and prove the long standing space lower bound of OðnÞ. The BiWFA edlib. In addition, for low error rates, memory-efficient WFA-based answers the pressing need for sequence alignment methods capable methods require a minimal amount of memory. to scaling to genome-scale alignments and full pangenomes. Most notably, BiWFA execution times are very similar, or even 3.3 Evaluation on real data better, than those of the original WFA (despite BiWFA requiring Figure 2 shows the performance results obtained for all the eval- 2954 and 607 less memory when aligning ultra-long MinION uated algorithms in terms of execution time and consumed memory. and PromethION sequences, respectively). This result can be better BiWFA uses many times less memory than other methods. In par- understood considering the memory inefficiencies that the original ticular, when aligning ultra-long ONT sequences (Fig. 2B), BiWFA WFA experiences when using a large memory footprint. As the requires between 68  93 less memory compared to wfalm and sequence’s length and error increases, the original WFA uses a sub- WFA low-memory modes. Furthermore, BiWFA uses 3.5 less stantially larger memory footprint, putting a significant pressure on memory compared to the efficient recursive mode from wfalm (most the memory hierarchy of the processor. Due to the pervasive mem- memory-efficient gap-affine algorithm to date). ory inefficiencies of modern processors executing memory intensive At the same time, BiWFA proves to be one of the fastest imple- applications, the original WFA’s performance is severely deterio- mentations aligning long sequences. Using ultra-long sequences, our rated when aligning long sequence datasets (like those from method is 25:7 faster than wfalm’s recursive mode. Moreover, Nanopore presented in the evaluation). In contrast, BiWFA relieves BiWFA’s execution times are similar to those of BitPal (sometimes this memory pressure using a minimal memory footprint. As a A Maximum memory consumption ONT PromethION reads vs CHM13 v1.1 > 10 kbps ONT Ultra Long > 500 kbps 1 TB 100 GB 10 GB 1 GB 100 MB 10 MB 1 MB B Execution time ONT PromethION reads vs CHM13 v1.1 > 10 kbps ONT Ultra Long > 500 kbps 1,000.0 Seconds 10.0 0.1 ) ) ) ) ce nly sse igh ed low alm low sive FA ce nly sse igh ed low alm low sive FA i s tan re o tz2− FA−h FA−m FA− wf alm− ecur BiW istan re o tz2− FA−h FA−m FA− wf alm− ecur BiW it d sc o ex W W W f w m−r it d sc o ex W W W wf m−r ed pal ( sw2− a l ed pal ( sw2− l ( lib bit k w f ( lib bit k wfa ed ed Algorithm Fig. 2. Experimental results from the execution of BiWFA and other state-of-the-art implementations aligning long sequences. Figure shows (A) memory consumption and (B) execution time per sequence aligned. A vertical line on each panel separates algorithms that use simpler penalty models or can only compute the alignment score (i.e. edlib and bitpal) from those that compute the full gap-affine alignment Optimal gap-affine alignment in O(s) space 7 result, BiWFA is able to balance out the additional work induced by Daily,J. et al. (2015) A work stealing based approach for enabling scalable op- BiWFA’s recursion, delivering a performance on-par with the origin- timal sequence homology detection. J. Parallel Distributed Comput., 79, al WFA. 132–142. We have presented the BiWFA using gap-affine penalties. Durbin,R. et al. (1998) Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press, Cambridge. Nevertheless, these very same ideas can be translated directly into Eizenga,J.M. and Paten,B. (2022) Improving the time and space complexity of other distances like edit, linear gap or piecewise gap-affine. the WFA algorithm and generalizing its scoring. bioRxiv. Moreover, it can be easily extended to semi-global alignment (a.k.a. Farrar,M. (2007) Striped Smith–Waterman speeds database searches six times ends-free, glocal, extension or overlapped alignment) by modifying over other SIMD implementations. Bioinformatics, 23, 156–161. the initial conditions and termination criterion. At the same time, Garrison,E. and Marth,G. (2012) Haplotype-based variant detection from the BiWFA retains the strengths of the original WFA: no restrictions short-read sequencing. arXiv preprint arXiv:1207.3907. on the sequences’ alphabet, preprocessing steps, nor prior estimation Gotoh,O. (1982) An improved algorithm for matching biological sequences. J. of the alignment error. Mol. Biol., 162, 705–708. Due to the simplicity of the WFA’s computational pattern, Jones,N. et al. (2004) An Introduction to Bioinformatics Algorithms. MIT BiWFA’s core functions can be easily vectorized to fully exploit the Press, Cambridge. capabilities of modern SIMD multicore processors. Our implemen- Kiełbasa,S. et al. (2011) Adaptive seeds tame genomic sequence comparison. Genome Res., 21, 487–493. tation, relies on the automatic vectorization capabilities of modern Koren,S. et al. (2017) CANU: scalable and accurate long-read assembly via compilers. As a result, the BiWFA implementation can exploit the adaptive k-mer weighting and repeat separation. Genome Res., 27, SIMD capabilities of any processor supported by modern compilers, 722–736. without rewriting any part of the source code. Li,H. (2013) Aligning sequence reads, clone sequences and assembly contigs Genomics and bioinformatics methods will continue to rely on with bwa-mem. arXiv preprint arXiv:1303.3997. sequence alignment as a core and critical component. BiWFA paves Li,H. (2018) Minimap2: pairwise alignment for nucleotide sequences. the way for the development of faster and more accurate tools that Bioinformatics, 34, 3094–3100. can scale with longer and noisy sequences using a minimal amount Loving,J. et al. (2014) BitPAl: a bit-parallel, general integer-scoring sequence of memory. In this way, we expect BiWFA to enable efficient se- alignment algorithm. Bioinformatics, 30, 3166–3173. quence alignment at genome-scale in years to come. Marco-Sola,S. et al. (2012) The gem mapper: fast, accurate and versatile align- ment by filtration. Nat. Methods, 9, 1185–1188. Marco-Sola,S. et al. (2021) Fast gap-affine pairwise alignment using the wave- front algorithm. Bioinformatics, 37, 456–463. Acknowledgements McKenna,A. et al. (2010) The genome analysis toolkit: a mapreduce frame- The authors thank Ragnar Groot Koerkamp and the anonymous reviewers work for analyzing next-generation DNA sequencing data. Genome Res., for making useful suggestions and contributing to improving the manuscript. 20, 1297–1303. Miga,K. H. and Wang,T. (2021) The need for a human pangenome reference sequence. Annu. Rev. Genomics Hum. Genet., 22, 81–102. Funding Myers,E. W. (1986) An O(ND) difference algorithm and its variations. Algorithmica, 1, 251–266. This research was supported by the European Union Regional Development Myers,E. W. and Miller,W. (1988) Optimal alignments in linear space. Fund within the framework of the ERDF Operational Program of Catalonia Bioinformatics, 4, 11–17. 2014-2020 with a grant of 50% of total cost eligible under the DRAC project Needleman,S. B. and Wunsch,C. D. (1970) A general method applicable to the [001-P-001723] and Lenovo-BSC Contract-Framework Contract (2022). It was search for similarities in the amino acid sequence of two proteins. J. Mol. also supported by the Ministerio de Ciencia e Innovacion MCIN AEI/10.13039/ Biol., 48, 443–453. Rahn,R. et al. (2018) Generic accelerated sequence alignment in seqan using 501100011033 and NextGenerationEU/PRTR under contracts PID2020- vectorization and multi-threading. Bioinformatics, 34, 3437–3445. 113614RB-C21, PID2019-107255GB-C21, and TED2021-132634A-I00, by Rodrıguez-Martın,B. et al. (2017) Chimpipe: accurate detection of fusion the Generalitat de Catalunya GenCat-DIUiE (GRR) [contracts 2017-SGR-313, genes and transcription-induced chimeras from RNA-seq data. BMC 2017-SGR-1328 and 2017-SGR-1414]. M.M. was partially supported by the Genomics, 18, 7–17. Spanish Ministry of Economy, Industry and Competitiveness under Ramon y Rognes,T. and Seeberg,E. (2000) Six-fold speed-up of Smith–Waterman se- Cajal fellowship number RYC-2016-21104. S.M.-S. was supported by Juan de quence database searches using parallel processing on common microproc- essors. Bioinformatics, 16, 699–706. la Cierva fellowship grant IJC2020-045916-I funded by MCIN/AEI/10.13039/ Simpson,J. et al. (2009) ABYSS: a parallel assembler for short read sequence 501100011033 and by ‘European Union NextGenerationEU/PRTR’. B.P. and data. Genome Res., 19, 1117–1123. J.M.E. were supported, in part, by the United States National Institutes of Smith,T. F. and Waterman,M. S. (1981) Comparison of biosequences. Adv. Health [award numbers: R01HG010485, U01HG010961, OT2OD026682, Appl. Math., 2, 482–489. OT3HL142481 and U24HG011853]. E.G. was supported by NIH/NIDA  Sosic,M. and Sikic,M. (2017) EDLIB: a C/Cþþ library for fast, exact sequence U01DA047638 and NSF PPoSS Award #2118709. A.G. acknowledges Dr. alignment using edit distance. Bioinformatics, 33, 1394–1395. Nicole Soranzo’s efforts to establish a pangenome research unit at the Human Suzuki,H. and Kasahara,M. (2017) Acceleration of nucleotide semi-global alignment with adaptive banded dynamic programming. bioRxiv. Technopole in Milan, Italy. Suzuki,H. and Kasahara,M. (2018) Introducing difference recurrence relations Conflict of Interest: none declared. for faster semi-global alignment of long sequences. BMC Bioinformatics, 19, 33–47. Wozniak,A. (1997) Using video-oriented instructions to speed up sequence comparison. Bioinformatics, 13, 145–150. References Xia,Z. et al. (2021) A review of parallel implementations for the Altschul,S. et al. (1990) Basic local alignment search tool. J. Mol. Biol., 215, Smith–Waterman algorithm. In: Interdisciplinary Sciences: Computational 403–410. Life Sciences, pp. 1–14. Bowden,R. et al. (2019) Sequencing of human genomes with nanopore tech- Zhang,Z. et al. (2000) A greedy algorithm for aligning DNA sequences. nology. Nat. Commun., 10, 1869. J. Comput. Biol., 7, 203–214. Daily,J. (2016) Parasail: SIMD C library for global, semi-global, and local Zhao,M. et al. (2013) SSW library: an SIMD Smith–Waterman C/Cþþ library pairwise sequence alignments. BMC Bioinf., 17, 1–11. for use in genomic applications. PLoS ONE, 8, e82138.