Bioinformatics, 37(4), 2021, 456–463 doi: 10.1093/bioinformatics/btaa777 Advance Access Publication Date: 11 September 2020 Original Paper Sequence analysis Fast gap-affine pairwise alignment using the wavefront algorithm 1,2, Santiago Marco-Sola *, Juan Carlos Moure2, Miquel Moreto1,3 and Antonio Espinosa2 1 Department of Computer Sciences, Barcelona Supercomputing Center, Barcelona 08034, Spain, 2Departament d’Arquitectura de Computadors i Sistemes Operatius, Universitat Autònoma de Barcelona, Barcelona 08193, Spain and 3Departament d’Arquitectura de Computadors, Universitat Politècnica de Catalunya, Barcelona 08034, Spain *To whom correspondence should be addressed. Associate Editor: Peter Robinson Received on April 8, 2020; revised on July 22, 2020; editorial decision on August 26, 2020; accepted on September 1, 2020 Abstract Motivation: Pairwise alignment of sequences is a fundamental method in modern molecular biology, implemented within multiple bioinformatics tools and libraries. Current advances in sequencing technologies press for the devel- opment of faster pairwise alignment algorithms that can scale with increasing read lengths and production yields. Results: In this article, we present the wavefront alignment algorithm (WFA), an exact gap-affine algorithm that takes advantage of homologous regions between the sequences to accelerate the alignment process. As opposed to trad- itional dynamic programming algorithms that run in quadratic time, the WFA runs in time O(ns), proportional to the read length n and the alignment score s, using Oðs 2 Þ memory. Furthermore, our algorithm exhibits simple data dependencies that can be easily vectorized, even by the automatic features of modern compilers, for different archi- tectures, without the need to adapt the code. We evaluate the performance of our algorithm, together with other state-of-the-art implementations. As a result, we demonstrate that the WFA runs 20–300 faster than other methods aligning short Illumina-like sequences, and 10–100 faster using long noisy reads like those produced by Oxford Nanopore Technologies. Availability and implementation: The WFA algorithm is implemented within the wavefront-aligner library, and it is publicly available at https://github.com/smarco/WFA. Contact: santiagomsola@gmail.com 1 Introduction development of faster algorithms that can keep up with the data pro- duction pace and scale with longer read lengths produced by newer Pairwise alignment of biological sequences is a core component of sequencing technologies. many bioinformatics tools. In genomics, it is an essential building In its most common formulation, the pairwise alignment prob- block for read mapping (Langmead and Salzberg, 2012; Li, 2013; lem is solved using some variation of the Needleman–Wunsch (NW) Marco-Sola et al., 2012), variant detection (DePristo et al., 2011), algorithm (Needleman and Wunsch, 1970) for gap-linear penalties de novo genome assembly (Simpson et al., 2009), multiple sequence or the Smith–Waterman–Gotoh (SWG) algorithm (Gotoh, 1982) for alignment (Notredame et al., 2000) and many other methods. gap-affine penalties. These solutions are based on dynamic program- Moreover, it can be adapted to different protocols and applications ming (DP) and consist on first computing the recurrence equations (e.g. DNA-seq, RNA-seq or bisulfite analysis) using different dis- on a DP matrix, and then tracing back the optimal alignment. As a tance functions (e.g. mismatch, edit, gap-linear or gap-affine), modi- result, they require quadratic time and memory on the length of the fying the alignment scope (e.g. global, semi-global or local sequences. alignment) or adjusting the scoring penalties. Nevertheless, all these Over the years, these algorithms have been extensively studied, variations share the same algorithmic core. and many different techniques and heuristics have been developed to Over the last decade, the volume of biological data generated by accelerate DP computations. The most successful optimizations sequencing technologies has increased exponentially. Besides, third- focus on vectorization strategies (Daily, 2016) (including bit-parallel generation technologies, such as those developed by PacBio and encodings Loving et al., 2014), banded computations (Suzuki and Oxford Nanopore Technologies (ONT), can yield sequences longer Kasahara, 2017) and cut-off techniques (Gusfield, 1997). than 20K bases; more than 100 longer than those produced by Notwithstanding these efforts, the quadratic execution time of clas- the extensively used Illumina sequencers. This presses for the sical approaches quickly becomes the bottleneck, and these methods C The Author(s) 2020. Published by Oxford University Press. V 456 This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Gap-affine wavefront alignment 457 fail to scale with longer read lengths. Besides, intrinsic dependencies 1.2 Our contribution on the DP recurrences limit the effectiveness of vectorization In this article, we present the wavefront alignment algorithm (WFA) approaches; what is more, they require ad hoc implementations for for exact pairwise alignment of sequences. This novel alignment each different SIMD instruction set. In many cases, practical imple- method progressively computes partial alignments of increasing mentations tend to apply heuristics that can miss the optimal align- score until the optimal solution is reached. In this way, the WFA al- ment at the expense of faster execution times. Interestingly, all these gorithm exploits similarities between sequences to accelerate the algorithms are oblivious to the similarities between the sequences. computation of the optimal alignment. As a result, its time complex- That is, they perform the same computations whether the sequences ity O(ns) depends on the sequence length n and the optimal align- are highly divergent or completely identical. ment score s. This allows scaling with sequence length provided that the error-rate between the sequences remains moderate. We demon- 1.1 Previous work strate that the WFA algorithm outperforms other state-of-the-art Over the years, alignment algorithms for the edit distance have been methods, while requiring less memory. largely studied and optimized due to its simplicity. As a result, many efficient libraries and tools like Edlib ( Sosic and  Sikic, 2017) [imple- menting Myers’s bit-vector algorithm (Myers, 1999)], or DAligner 2 Materials (Myers, 2014) [implementing Myers’s O(ND) difference algorithm In the following sections, we describe the main concepts behind the (Myers, 1986)] are now available. However, the edit distance fails WFA algorithm. First, we introduce the concept of furthest-reaching to capture critical properties of the difference between biological points and wavefronts, and then we redefine the recurrent equations sequences. For that reason, gap-affine distance is commonly pre- ferred, and thus many research efforts have been invested on acceler- of the SWG algorithm in terms of wavefronts. Finally, we prove that ating the SWG algorithm. the first wavefront that reaches the bottom-right corner of the DP Many of the SWG optimization techniques are focused on matrix corresponds to the optimal pairwise alignment between the exploiting intra-sequence parallelism using tailored vectorization sequences. strategies. In this way, Wozniak (1997) first proposed to compute DP cells along the minor diagonal in parallel. Then, Rognes (Rognes 2.1 Pairwise global alignment of sequences and Seeberg, 2000) proposed to vectorize along the query sequence, Let the query q ¼ q0 q1 . . . qn1 and the text t ¼ t0 t1 . . . tm1 be introducing the query profile technique. Later on, Farrar (2007) pro- strings of length jqj ¼ n and jtj ¼ m, respectively (Table 1). We de- posed to reorganize computations in a striped fashion. Farrar’s fine the pairwise global alignment problem in terms of computing method dominated over all other approaches until Daily presented the alignment [or edit path, Myers (1986)] from (0, 0) to (n, m) with the prefix-scan method (Daily et al., 2015), which in practice out- minimum penalty score, allowing matching bases, substitutions and performs all previous methods. Meanwhile, focusing on further gaps (i.e. insertions and deletions). Under the gap-affine model, the exploiting instruction-level parallelism, the BitPAl algorithm alignment penalty score is computed based on fa; x; o; eg, where a (Loving et al., 2014) was presented as a natural extension of and x correspond to the penalty of matching or mismatching two Myers’s bit-vector algorithm restricted to gap-linear penalties. bases, respectively, and the gap-penalty function is expressed as the Although it requires generating tailored code for each different set of penalties, it represents an interesting method for computing the linear function gðnÞ ¼ o þ n  e (where n is the length of the gap). As alignment distance (i.e. not the full alignment) if all the penalties are stated before, this problem is commonly solved using some variation fixed and known in advance. of the SWG algorithm (Gotoh, 1982). Eq. 1 shows the recurrence Recently, Suzuki and Kasahara (2018) proposed a differential relations of the DP matrix used in the SWG algorithm. encoding of the DP matrix that effectively reduces the size required 8 < Iv;h ¼ minfMv;h1 þ o þ e; Iv;h1 þ eg to store each DP cell to further exploit SIMD operations. Combined Dv;h ¼ minfMv1;h þ o þ e; Dv1;h þ eg : (1) with an adaptive-band, that adjusts its size according to the scores : Mv;h ¼ minfIv;h ; Dv;h ; Mv1;h1 þ sðqv1 ; th1 Þg computed so far, they developed the Gaba library (Suzuki and Kasahara, 2017). Because of its good performance, this algorithm Despite its general formulation, in practice, the majority of was re-implemented in the KSW2 library, which is currently used applications use strictly positive penalties, i.e. x; o; e > 0. And, as in within the core of the Minimap2 mapper (Li, 2018). other methods, we focus on the match-mismatch model where the Over the past decade, many libraries implementing efficient mismatch between every two different bases always has the same SWG algorithms have been developed due to its critical role in sev- penalty, x. That is, function s(v, h), for 0  v < n and eral bioinformatics tools and pipelines. The SSW library (Zhao 0  h < m, evaluates to a if qv ¼¼ th and x otherwise. In addition, et al., 2013) was one of the first libraries to offer an efficient imple- our method fixes the match score to zero (a ¼ 0) to enable exploiting mentation of the SWG using Farrar’s striped algorithm. Later on, similarities between the sequences and accelerate the alignment pro- the Parasail library (Daily, 2016), that implements several intra- cess. As we demonstrate later, this allows defining a much simpler sequence vectorization strategies, was published. But also, the SeqAn (Rahn et al., 2018) developers have put many efforts into and faster method for gap-affine alignment. incorporating high-performance and parallel implementations of the SWG among its many algorithms. At the same time, all these libra- Table 1. Notations ries implement additional optimizations that can further enhance the presented methods. For instance, the integer saturation strategy, Symbol Description which first computes the alignment using a small integer type, and recomputes the alignment in case an overflow/underflow is detected. q, t Strings: q ¼ q0 q1 . . . qn1 ; t ¼ t0 t1 . . . tm1 In addition, these libraries implement several coarse-grain parallel- p ¼ fa; x; o; eg Gap-affine penalty scores ization techniques that allow performing various independent align- M, I, D SWG alignment matrices ments in parallel. But also, they include tailored implementations F s;k Furthest-reaching (f.r.) point on diagonal k with score s for specific SIMD instruction sets than can yield shorter execution M ~ s;k Offset in the diagonal to the corresponding f.r. point F s;k times. I~s ; D~s ; M ~s Components of the wavefront WFs Last, it should be noted that there are many other lines of research WFs Wavefront of score s for accelerating the SWG involving hardware accelerators; including ~ hi M Index of the rightmost diagonal in the component Graphics Processing Units (Chacón et al., 2014; Liu et al., 2013), ~ lo M Index of the lowest diagonal in the component Field-Programmable Gate Arrays (Li et al., 2007), Cell Broadband WFmin Minimum wavefront length Engine (Szalkowski et al., 2008), Xeon Phi processor (Liu et al., 2014) WFdiff Maximum difference distance and even custom hardware designs (Hasan et al., 2008). 458 S.Marco-Sola et al. 2.2 Wavefront furthest-reaching diagonals 2.3 Wavefront algorithm Restricted to the edit distance, the so-called diagonal-transition Using Lemma 2.1 and staring with M ~ 0;0 ¼ 0, the WFA algorithm algorithms (Landau and Vishkin, 1989; Ukkonen, 1985) exploited (Algorithm 1) progressively computes wavefronts of increasing score the fact that DP diagonals have monotonically increasing scores. In until the cell (n, m) is reached. First, for each score s, the algorithm particular, Myers’ O(nd) algorithm (Myers, 1986), implemented at the core of the Linux diff-tool, computes the longest common subse- quence of two sequences taking advantage of stretches of matching characters along the diagonals. Our algorithm extends these ideas to the gap-affine model. For each score s and diagonal k, the furthest-reaching (f.r.) point F s;k denotes the DP-cell, on diagonal k and with score s, that is more far-off from the beginning of the diagonal. Note that a point p ¼ ðv; hÞ, in the diagonal k ¼ h  v, is further than other p0 ¼ ðv  1; h  1Þ (p0 < p) in the same diagonal. Then, we define ~ s;k and M I~s;k ; D ~ s;k as the offset in the diagonal to the f.r. point F s;k for each of the three SWG matrices I, D and M, respectively. So, for a given score s, we define the s-wavefront (WFs) as the set of all the f.r. points with score s; that is, the set of offsets I~s;k ; D ~ s;k and M ~ s;k , for all k. Likewise, we denote I~s ; D~s and M ~ s , the components of the wavefront WFs. Assuming each component of a wavefront is repre- sented using a vector of offsets centred over the main diagonal (k ¼ 0), let M ~ hi denote the index of the rightmost diagonal in the component, and M ~ lo the index of the lowest. Our goal is to compute the minimum s such that any of the f.r. points of WFs reaches (n, m). We notice that for any s, the f.r. points of WFs can only be originated from points whose score is s – o, s – e or s – x; or from a previous point with score s solely followed by matches along the diagonal. Considering only insertions, deletions and mis- matches; we can redefine Eq. 1 in terms of offsets to f.r. points (Eq. 2). ( ) M ~ soe;k1 ðOpen insertionÞ I~s;k ¼ max þ1 I~se;k1 ðExtend insertionÞ ~ s;k following matching characters along the ( ) extends the points M ~ diagonals, using WF_EXTEND(M ~ s ; q; t) (Algorithm 2). Then, it D ~ s;k ¼ max M soe;kþ1 ðOpen deletionÞ ~ se;kþ1 checks whether any of the resulting f.r. points of wavefront WFs D ðExtend deletionÞ (2) 8 9 reaches (n, m). If not, the algorithm proceeds to compute the next > ~ > wavefront WFsþ1 using WF_NEXT(D; ~ I;~ M; ~ q; t; s) (Algorithm 3, < M sx;k þ 1 ðSubstitutionÞ > > = which applies Eq. 2), and iterates again.In the case of global align- M ~ s;k ¼ max I~ ðInsertionÞ ; > s;k > ment, the algorithm starts with a unitary wavefront (M ~ 0;0 ). As the > :D ~ > ; s;k ðDeletionÞ with initial condition M ~ 0;0 ¼ 0 Essentially, Eq. 2 states that we can compute the set of f.r. points of WFs using WFso ; WFse and WFsx . For example, in order to compute the f.r. point D ~ s;k , we only have to consider which f.r. point brings the resulting offset further: a newly opened deletion M ~ soe;kþ1 , or an extended one D ~ se;kþ1 (both coming from one diagonal above). Similarly, as to compute the f.r. point M~ s;k , we only need to compare offsets between the f.r. points cor- responding to an insertion I~s;k , a deletion D ~ s;k and a mismatch from M ~ sx;k . Then, we consider matches, and we take advantage of the fact that a ¼ 0. For that, we extend all the previously com- puted points by Eq. 2, as far as possible, following matching char- acters along the diagonal. LEMMA 2.1. The resulting points of applying Eq. 2 and extending them, following matching characters along the diagonal, are the f.r. points for algorithm iterates, the length of the wavefronts increases according the score s, that is WFs. to Eq. 3. That is, each new wavefront grows to span over one more diagonal on each end (i.e. M~ hi and M ~ lo ) compared to the wavefronts Proof. Suppose that point p 2 WFs (originated by point q from a previ- s s it depends on. As a result, the size of each subsequent wavefront ous WFs0 ; s0 ¼ s  p; p 2 fx; o; eg) is not the f.r. point on the diagonal k. increases proportional to the alignment score between the sequen- Then, there must exist a point p0 that is f.r. on that diagonal ðp0 > pÞ. ces. Hence, the algorithm requires Oðs2 Þ memory to store all For that p0 , there must exist a point q0 (with a score s0 ) that generated p0 wavefronts. (perhaps, after extension by accounting matches on the diagonal). n o However, q0 cannot be the f.r. point on the diagonal (i.e. does not belong ~ hi ¼ I~hi ¼ D ~ hi ¼ max M~ hi ; M ~ hi ~hi ; D~ hi þ 1 M ; I to WFs0 ), otherwise, it would have generated p. Therefore, q0 < q and s s s sx soe se se n o (3) thus, it follows that p0  p, in the best case (f.r point), p0 ¼ p. h ~ lo ~lo ~ lo ~ lo ~ lo ~se ~ se M s ¼ I s ¼ D s ¼ max M sx M soe I s ; D s  1: It follows from Lemma 2.1 that the optimal alignment corresponds to the sequence of wavefronts from WF0 to WFs, being WFs the wavefront Also, note that extending a wavefront (WF_EXTEND function) with the smallest s that contains a f.r. point that reaches (n, m). is bounded by the number of diagonal matching characters Gap-affine wavefront alignment 459 (maxfn; mg) and the length of the wavefront. Similarly, the function WF_NEXT computes each next wavefront in time proportional to the wavefront length. Therefore, the running time of the WFA algo- rithm, to compute an alignment of score s, is bounded in the worst case by Oðmaxfn; mg  sÞ, or O(ns) assuming that the sequences have the same length. Figure 1 shows an example of computing a global gap-affine alignment using a DP matrix and the WFA algorithm. The figure demonstrates that the WFA only needs to compute a minimal num- ber of DP cells to find the optimal alignment. The WFA explores the cells of the DP matrix by increasing score, like seeking for the min- imum cost path from (0, 0) to (n, m). In this way, the algorithm ef- fectively computes an adaptive band without previous knowledge of the alignment error between the sequences. Moreover, it performs this seemingly irregular computation using simple and regular opera- Fig. 1. (A) DP matrices (i.e. M, D and I) depicting the only cells computed by the WFA in order to align q ¼ “GATACA” against t ¼ “GAGATA” using the penalties tions over the wavefronts. Furthermore, the algorithm only requires fx ¼ 4; o ¼ 6; e ¼ 2g. The arrows show the f.r. points of the main diagonal in ma- to store a few wavefronts before reaching the bottom right cell of the ~ 0;0 ; M trix M (i.e. M ~ 4;0 and M ~ 8;0 ). (B) Equivalent representation using wavefronts DP matrix. Considering the penalties used in Figure 1, the only score (i.e. WF0, WF4 and WF8) values that can appear in the DP matrix, before reaching the optimal alignment at s ¼ 8, are s ¼ 0 and s ¼ 4. For that reason, the WFA only needs to store the wavefronts WF0, WF4 and WF8. And, unlike the DP matrix, these wavefronts avoid storing runs of cells with the same score (byproduct of the matches along the diagonal). Once the algorithm computes a wavefront that reaches (n, m), and therefore, the optimal alignment between the sequences, it can retrieve back the path that leads from (0, 0) to (n, m) (i.e. back- trace). As opposed to DP-based algorithms, the WFA’s backtrace is performed across the wavefronts’ offsets instead of using the DP ma- trix scores. On each step of the backtrace, the function determines which f.r. point, from the previous wavefronts, originated the cur- rent offset. The difference between the actual offset and the source is the total amount of matching characters between the two positions. For example, in Figure 1, offset M ~ 8;0 ¼ 6 has been generated from ~ 4;0 ¼ 4. The difference between the two offsets corre- the f.r. point M sponds to the matches from (5, 5) to (6, 6), both included. 2.4 Adaptive wavefront reduction As the WFA algorithm iterates, wavefronts of increasing score span over more and more diagonals. While some of the f.r. points quickly distance for all the f.r. points of the wavefront carries a non- negligible toll on the total running time of the algorithm. For that reason, we define a minimum wavefront length WFmin required to trigger the reduction (WF_REDUCE). In practice, the WFA-Adapt heuristic successfully removes dis- tant f.r. points unlikely to be part of the optimal alignment. advance towards the solution, others are left behind on their diag- Therefore, the span of the wavefront is focused on the diagonals onal as they are unlikely to lead to the optimal solution. As a conse- closer to the alignment solution, further reducing the amount of quence, the WFA algorithm invests a substantial amount of time memory required by the algorithm. In a way, similar pruning techni- processing unpromising paths. For that reason, we propose a heuris- ques are implemented into DP-based algorithms. These banded tic version of the algorithm (WFA-Adapt) that removes f.r. points methods avoid computing paths beyond a band of diagonals. Unlike from outer diagonals that are extremely far behind compared to those methods, our wavefront reduction method prunes based on other points in the same wavefront. the potential of the diagonal to lead to the optimal solution, without In this way, the WFA-Adapt performs a wavefront pruning previous knowledge of the error between the sequences whatsoever. (Algorithm 4) after performing each wavefront extension. It deter- Nonetheless, it is important to highlight that the WFA-Adapt mines the distance dk of each diagonal f.r. point to the (n, m) cell. method is no longer exact and might miss some optimal solutions; Then it computes the minimum dmin ¼ mink fdk g and discards outer unlike the WFA. In practice, for a reasonable selection of diagonal points that are more than WFdiff cells behind it (i.e. ðWFmin ; WFdiff Þ, this is rarely the case, even when aligning long dk  dmin > WFdiff ). Despite its simplicity, computing each dk noisy sequences. 460 S.Marco-Sola et al. 2.5 Efficient computation of wavefronts 3 Results As opposed to traditional DP-based algorithms, the core functions We evaluate the performance of our implementation of the WFA of the wavefront algorithm depict simple data dependencies (Fig. 2). and WFA-Adapt algorithms, together with other widely used libra- This allows the implementation of effective fine-grain parallelization ries that implement optimized versions of well-known pairwise techniques to accelerate the computation of each wavefront. alignment algorithms. In the case of the WF_NEXT function, the automatic vectoriza- tion features of modern compilers can easily detect the dependencies between wavefronts and transparently issue SIMD instructions to 3.1 Experimental setup significantly accelerate the computation. This represents a major ad- We select the widely used Parasail library (Daily, 2016) that imple- vantage over other methods that rely on specific SIMD intrinsics ments the diagonal method by Wozniak (Wozniak, 1997) and custom implementations. These tools not only have limited (ParaDiag), the striped approach proposed by Farrar (Farrar, 2007) portability to different computer architectures but also require com- (ParaStrip), the prefix-scan algorithm (Daily et al., 2015) (ParaScan) plicated, time-consuming and error-prone implementations. In con- and the classical banded heuristic method (ParaBand). Also, we evaluate the gap-affine implementation of the SeqAn library (Rahn trast, our approach allows transparent vectorization of the et al., 2018), the bit-parallel implementation for non-unitary penal- algorithm for any SIMD instruction-set and architecture supported ties (Loving et al., 2014) (BitPAl) and the adaptive band algorithm by the compiler. from Suzuki (Suzuki and Kasahara, 2018) (Gaba). In addition, we In the same way, the function WF_EXTEND can be accelerated evaluate the performance of the highly optimized SIMD algorithms by means of exploiting bit-level parallelism techniques. That is, in- KSW2-Z2 (ksw2_extz2_sse), and KSW2-D2 (ksw2_extd2_sse), stead of comparing characters along the diagonal one at a time, we from the KSW2 library, used within Minimap2 (Li, 2018). On top perform comparisons in blocks of eight characters (i.e. packed in a of that, we present the performance results of the Edlib ( Sosic and 64-bit computer word). In this way, the algorithm performs parallel  Sikic, 2017) library as a base line for comparison. Because Edlib is comparisons and extends each diagonal up to eight positions per restricted to compute the edit-distance alignment, its execution is re- step. Taking into account that spurious matches between blocks of markably fast. eight characters are very unlikely to occur, this technique prevents For the purpose of the evaluation, all the methods have been con- the inner loop of WF_EXTEND from iterating more than once in figured to generate global alignments. These methods are grouped in most cases. Otherwise, it means that the function has detected a three categories: ‘Gap-affine (Exact)’ (WFA, KSW2-Z2, SeqAn, stretch of matches which, in turn, favors the overall performance of ParaStrip, ParaScan and ParaDiag) for exact algorithms that use the algorithm. gap-affine penalties, ‘Gap-affine (Banded)’ (WFA-Adap, Gaba and Note that, depending on the penalties selected, some (or even all) ParaBand) for approximate algorithms that use gap-affine penalties, of the wavefronts required to compute a wavefront might not exist. and ‘Others’ (Edlib, BitPAl and KSW2-D2) for exact methods that That means that there is no possible combination of penalties that use other penalty models. Note that the gap-affine banded methods can generate the current score. In those cases, the number of opera- are approximated and not guaranteed to find the optimal alignment. tions involved in WF_NEXT gets reduced. Hence, our implementa- Similarly, note that the algorithms ParaBand and BitPAl only com- tion contains different specializations of this function to accelerate pute the alignment score (not the full alignment) and the latter is those cases where some source wavefronts are not needed, and thus restricted to linear-gap penalties. the complexity can be reduced. We discarded other methods from the evaluation as their run- Notably, the WFA encodes the offsets to the f.r. points on each ning time was exceedingly long or because their recall was substan- diagonal, not the actual score values. Therefore, the representation tially below par. This is the case of other KSW2 algorithms (like the range of each wavefront element is bounded by the length of the approximate versions), other edit distance methods [e.g. BPM sequences, and not by the maximum score possible. In most cases, (Myers, 1999) or DAligner (Myers, 2014)] or custom implementa- the WFA can conveniently encode offsets using 16-bit integers (for tions for hardware accelerators (which fall beyond the scope of this sequences up to 216 characters), or even 8-bit integers (for Illumina- article). like sequences, no longer than 255 bases). This succinct representa- We choose two real datasets: 100 000 sequences from an tion not only allows saving memory but also allows using wider Illumina HiSeq 2000 of 100 bp (Accession number ERX069505), and 25 000 ONT MinION sequences from (Bowden et al., 2019) SIMD instructions that can operate over more elements at once. (Accession numbers ERR3278877–ERR3278886). Also, we simu- In turn, this further enhances the efficiency of our WFA late several datasets of various lengths (i.e. 100, 1K, 10K and 100K implementation. bases), with different error rates (i.e. d ¼ 1, 2, 5, 10 and 20% error rate) by means of randomly adding mismatches and indels. For every dataset, we select as many sequences as needed, so that each of them contains a total of 10 million bp. All tests are executed on an Intel Xeon Platinum 8160 equipped with 96 GB of RAM, running SuSE Linux Enterprise Server. For each method, we measure execution time and memory consumed. In addition, we verify the results using a basic SWG implementation in order to compute the recall of each method. 3.2 Evaluation on real data Table 2 shows the time performance results obtained for all the algorithms evaluated using both real and simulated datasets. On real datasets, the WFA performs many times faster than other methods. In particular, aligning HiSeq sequences, our method is 200–300 faster compared to traditional DP algorithms, and 20–40 faster than the adaptive-band methods. Similarly, when aligning ONT sequences, the WFA performs 28–200 and 6–7 times faster than DP algorithms and adaptive-band methods, respectively. Moreover, it is several times faster than methods that only compute the align- Fig. 2. (A) Dependencies between DP cells when computing another cell using gap- ment score. Furthermore, the WFA-Adapt refinement attains an affine penalties. (B) Dependencies between wavefronts as to compute one element of extra 1.6 speedup over the original WFA when aligning long ONT the next wavefront (Algorithm 3) reads. Gap-affine wavefront alignment 461 Table 2. Time performance of pairwise alignment algorithms 100 1K 10K 100K HiSeq ONT d ¼ 1% d ¼ 5% d ¼ 20% d ¼ 1% d ¼ 5% d ¼ 20% d ¼ 1% d ¼ 5% d ¼ 20% d ¼ 1% d ¼ 5% d ¼ 20% Gap-affine WFA 0.06 2.46 0.09 0.37 1.55 0.14 0.93 6.93 0.43 7.28 66.00 8.49 102.00 2542.00 (Exact) KSW2-Z2 2.42 15.23 2.41 2.79 3.06 16.43 16.47 16.64 188.40 188.40 189.00 2146.00 2136.00 2139.00 SeqAn 18.18 180.60 18.28 19.17 22.64 184.80 190.80 213.60 n/a n/a n/a n/a n/a n/a ParaStrip 77.40 530.40 77.40 75.60 70.80 603.60 602.40 572.40 n/a n/a n/a n/a n/a n/a ParaScan 12.03 69.00 11.85 11.89 11.99 75.00 75.00 75.00 746.40 747.00 747.00 n/a n/a n/a ParaDiag 17.80 130.20 17.53 17.58 17.66 141.60 141.60 141.60 n/a n/a n/a n/a n/a n/a Gap-affine WFA-Adap 0.07 1.51 0.09 0.50 2.20 0.16 0.68 2.75 0.17 0.68 3.00 0.26 1.53 6.97 (Banded) Gaba 1.26 error 1.28 1.33 error 0.68 0.74 error 0.60 0.81 error 0.72 0.76 error ParaBand 3.71 8.26 3.75 3.89 4.17 4.04 4.51 5.34 4.87 6.39 8.87 8.11 11.81 20.43 Others Edlib 1.75 2.28 1.73 1.82 1.94 1.95 2.03 2.73 3.67 4.50 7.39 6.95 19.34 error BitPAl 0.47 2.18 0.47 0.47 0.47 2.27 2.31 2.33 21.44 21.44 21.47 212.40 212.40 212.40 KSW2-D2 2.87 18.13 2.86 3.23 3.49 19.56 19.64 19.81 219.60 219.60 220.20 2385.00 2380.00 2452.00 Note: Running time in seconds obtained by the different pairwise alignment algorithms on real and simulated datasets. Penalties used for all tests fx ¼ 4; o ¼ 6; e ¼ 2g correspond to the ones used by BWA-MEM. Alignment algorithms have been grouped in three categories: ‘Gap-affine (Exact)’ for exact algorithms that use gap-affine penalties, ‘Gap-affine (Banded)’ for approximate algorithms that use gap-affine penalties, and ‘Others’ for exact algorithms that use other distance models (i.e. edit distance or linear-gap penalties). On simulated datasets, results are broken down by error rate. Experiments whose running time exceeded 2 hours are marked as ‘n/a’. Likewise, runs that consistently failed are labelled with ‘error’. Note that WFA-Adapt was executed using WFmin ¼ 10 and WFdiff ¼ 50. 3.3 Evaluation on simulated data turn, this is subject to the penalty scores chosen for the alignment. When aligning simulated datasets, WFA methods consistently out- Table 3 reflects how the performance of the WFA methods is perform other algorithms. Overall, WFA and adaptive-band meth- affected by different choices of penalty scores. Although the running ods are faster than classical DP algorithms. In particular, as the read times remain remarkably low, selecting penalties that increase the length increases, Parasail and SeqAn’s implementations require im- final alignment score can negatively impact the performance of the practical running times. Notably, BitPAl and ParaBand running algorithm; in this case, increasing the alignment time up to 2 in times remain reasonably low, but alas, their results are limited to the the worst case. alignment score. Altogether, these methods prove to be completely insensitive to the error between the sequences. For that reason, WFA 3.4 Memory footprint and WFA-Adapt methods achieve remarkable speedups when the Table 4 presents the overall memory consumption of the algorithms error rate remains moderate (i.e. d ¼ 1–5%). For instance, compared evaluated. These results experimentally confirm the reduced mem- to ParaBand, WFA methods run up to 13 and 78 faster for error ory footprint required by the WFA methods. Moreover, as the read rates of 1 and 5%, respectively; and up to 138 and 814 faster length and error rate increases, the WFA-Adapt method uses 19– compared to BitPAl for the same error rates. 36 times less memory than other methods. This not only results in On the other hand, WFA methods prove to be generally faster a significant reduction on the overall memory used but also relieves and more accurate than adaptive-band methods. Although these the pressure put in the memory hierarchy by traditional DP-based adaptive algorithms prove to be superior to classical approaches, algorithms. WFA and WFA-Adapt methods scale better with both the sequence length and the error rate. In particular, executions of the KSW2 li- brary consistently take between 2 and 510 longer than the WFA. 4 Discussion In the demanding scenario of aligning reads of 100K bases at just 1% of divergence rate, WFA-Adapt runs more than three orders of In this article, we presented a novel algorithm for gap-affine pair- magnitude faster than the KSW2 algorithms. It is important to note wise alignment. Our method exploits the similarities between the that KSW2-D2 uses 2-piece affine model (Gotoh, 1990), that is why sequences to deliver exact alignments while outperforming in time it is slower than KSW2-Z2 which uses the 1-piece/standard model. and memory other state-of-the-art algorithms. To this end, our algo- In general, as the error rate increases, adaptive-band methods ei- rithm computes alignments of increasing score using a succinct ther take longer to finish, like KSW2-Z2 and KSW2-D2, or become diagonal-transition representation. As a result, the WFA algorithm insensitive reporting suboptimal alignments, like Gaba. Table 3 runs in O(ns) time and memory; being n the sequence length, and s presents a summary of time and recall results obtained using two the alignment score between the sequences. broadly used penalty schemes; that is, those used by the ubiquitous To date, adaptive-band algorithms dominated over other mappers BWA-MEM (Li, 2013) and Bowtie2 (Langmead and approaches. In practice, these methods seek to heuristically delimit Salzberg, 2012). In the case of the Gaba algorithm, executions times the alignment path to avoid computing DP cells whose score is too seem to decrease as the sequence length increases. Nonetheless, this high. Compared to them, the WFA algorithm naturally computes is the result of aggressive heuristics that, in the end, cause a signifi- cells of the DP matrix by increasing score without introducing fur- cant drop in the total recall (5%) or fail to align noisy sequences ther complexities. Due to its simplicity, the WFA algorithm can be (e.g. aligning ONT sequences). Notably, only the WFA methods easily vectorized using SIMD instructions, as opposed to traditional could successfully align all datasets of 100K bases, whereas other DP-based algorithms. Moreover, data dependencies can be automat- methods took an unreasonable amount of time or failed. Moreover, ically understood by modern compilers in order to transparently WFA-Adapt could finish all the executions in less than 7 seconds issue SIMD instructions for any supported vectorial architecture. achieving the maximum possible recall; that is, reporting 100% of Furthermore, when aligning moderately long sequences (i.e. less the optimal alignments. than 255 bases), the WFA can encode diagonal offsets using 8-bit In practice, the WFA algorithm seems barely sensitive to the se- integers, which not only enhances SIMD performance but also fur- quence length, but to the alignment score between the sequences. In ther reduces the memory footprint. 462 S.Marco-Sola et al. Table 3. Accuracy of pairwise alignment algorithms HiSeq ONT l ¼ 1K, d ¼ 5% l ¼ 10K, d ¼ 5% Time Recall Time Recall Time Recall Time Recall BWA-MEM WFA 0.06 100 2.46 100 0.93 100 7.28 100 WFA-Adap 0.07 100 1.51 99.8 0.68 100 0.68 100 Gaba 1.26 100 error error 0.74 95.2 0.81 94.8 KSW2-Z2 2.42 100 15.23 100 16.47 100 188.4 100 KSW2-D2 2.87 100 18.13 100 19.64 100 219.6 100 Bowtie2 WFA 0.07 100 4.80 100 1.60 100 12.05 100 WFA-Adap 0.08 100 3.18 99.8 1.32 100 1.25 100 Gaba 1.26 100 error error 0.74 95.4 0.66 95.4 KSW2-Z2 2.40 100 15.21 100 16.46 100 190.8 100 KSW2-D2 2.91 100 18.16 100 19.66 100 220.8 100 Note: Time (in seconds) and recall (as percentage of the total sequences aligned) for error-sensitive algorithms. Two different penalty schemes were used: first corresponding to the default values used by BWA-MEM (i.e. x ¼ 4; o ¼ 6; e ¼ 2), and second matching Bowtie2’s end-to-end alignment scores (i.e. x ¼ 6; o ¼ 5; e ¼ 3). Table 4. Memory consumption of pairwise alignment algorithms l ¼ 100 l ¼ 1K l ¼ 10K HiSeq ONT d ¼ 1% d ¼ 20% d ¼ 1% d ¼ 20% d ¼ 1% d ¼ 20% WFA 0.5 3.1 0.4 0.4 1.4 4.2 6.2 269.6 KSW2-Z2 0.4 5 0.3 0.3 4.3 5 384.7 387.1 SeqAn 0.3 2.6 0.4 0.4 2.9 1.4 n/a n/a ParaStrip 6.7 7.9 68.5 15.6 6 13.1 n/a n/a ParaScan 6.8 7.9 68.5 15.4 6.2 12.8 196.1 202.6 ParaDiag 6.8 7.2 15.6 15.3 5.1 12 n/a n/a WFA-Adap 0.4 2 0.4 0.5 0.9 1.5 6.1 10.7 Gaba 0.7 Error 0.7 Error 4.3 Error 640.9 Error ParaBand 0.4 0.3 0.5 0.4 0.4 0.5 0.5 0.5 Edlib 0.3 0.4 0.3 0.3 0.3 0.3 0.4 0.4 BitPAl 0.4 0.5 0.5 0.3 0.3 0.4 0.6 0.7 KSW2-D2 0.3 4.7 0.4 0.3 4.3 4.6 384.8 388.5 Note: Memory usage (MB) of each alignment algorithm for different sequence length and error rate. In addition, the WFA algorithm poses no constrain on the alpha- Funding bet size, nor requires any preprocessing step in advance. Besides, the This research was supported by the European Unions’s Horizon 2020 progressive computation of alignments of increasing score allows Framework Programme under the DeepHealth project [825111], by the the WFA to implement exact cut-offs techniques. And yet, this method does not require any prior estimation of the alignment score European Union Regional Development Fund within the framework of the between the sequences whatsoever. Also, note that the optimizations ERDF Operational Program of Catalonia 2014-2020 with a grant of 50% of presented in this article are focused on intra-sequence paralleliza- total cost eligible under the DRAC project [001-P-001723], by the MINECO- tion. Consequently, all the experiments were executed using a single Spain (contracts TIN2017-84553-C2-1-R and TIN2015-65316-P) and by the thread. Nevertheless, the WFA can be used together with many Catalan government (contracts 2017-SGR-313, 2017-SGR-1328 and 2017- inter-sequence paralellization techniques in order to exploit the SGR-1414). M.M. was partially supported by the Spanish Ministry of multi-threading capabilities of modern processors. And, although Economy, Industry and Competitiveness under Ramón y Cajal fellowship the WFA algorithm has been presented using the gap-affine scoring number RYC-2016-21104. model, it can be easily adapted to more other scoring models, like the linear-gap model or the piece-wise affine model (Gotoh, 1990; Conflict of Interest: none declared. Miller and Myers, 1988). Likewise, it can be adapted to semi-global alignment by adjusting the initial conditions of Eq. 2 and the exit condition on Algorithm 1. But also, it can even be used for finding References overlaps and local alignments in the spirit of mappers like DAligner Bowden,R. et al. (2019) Sequencing of human genomes with nanopore tech- (Myers, 2014). nology. Nat. Commun., 10, 1–9. Undoubtedly, pairwise alignment will remain as a central and Chacón,A. et al. (2014) Thread-cooperative, bit-parallel computation of critical building-block of many bioinformatics applications. For that levenshtein distance on GPU. In Proceedings of the 28th ACM international reason, our algorithm represents a fast and scalable solution for conference on Supercomputing(ICS ’14). Association for Computing many sequence analysis tools to cope with the ever-increasing yields Machinery, New York, NY, USA, pp. 103–112. 10.1145/2597652. of sequencing technologies. In this way, the WFA approach will 2597677 pave the way for the design of better alignment algorithms in years Daily,J. (2016) Parasail: SIMD C library for global, semi-global, and local to come. pairwise sequence alignments. BMC Bioinformatics, 17, 81. Gap-affine wavefront alignment 463 Daily,J. et al. (2015) A work stealing based approach for enabling scalable op- Miller,W. and Myers,E.W. (1988) Sequence comparison with concave weight- timal sequence homology detection. J. Parallel Distrib. Comput., 79-80, ing functions. Bull. Math. Biol., 50, 97–120. 132–142. Myers,E.W. (1986) Ano (nd) difference algorithm and its variations. DePristo,M.A. et al. (2011) A framework for variation discovery and genotyp- Algorithmica, 1, 251–266. ing using next-generation DNA sequencing data. Nat. Genet., 43, 491–498. Myers,G. (1999) A fast bit-vector algorithm for approximate string matching Farrar,M. (2007) Striped smith–waterman speeds database searches six times based on dynamic programming. J. ACM (JACM), 46, 395–415. over other SIMD implementations. Bioinformatics, 23, 156–161. Myers,G. (2014) Efficient local alignment discovery amongst noisy long reads. Gotoh,O. (1982) An improved algorithm for matching biological sequences. In International Workshop on Algorithms in Bioinformatics. Springer, pp. J. Mol. Biol., 162, 705–708. 52–67. Gotoh,O. (1990) Optimal sequence alignment allowing for long gaps. Bull. Needleman,S.B. and Wunsch,C.D. (1970) A general method applicable to the Math. Biol., 52, 359–373. search for similarities in the amino acid sequence of two proteins. J. Mol. Gusfield,D. (1997) Algorithms on stings, trees, and sequences: computer sci- Biol., 48, 443–453. ence and computational biology. ACM Sigact News, 28, 41–60. Notredame,C. et al. (2000) T-coffee: a novel method for fast and accurate Hasan,L. et al. (2008) Hardware implementation of the smith-waterman multiple sequence alignment. J. Mol. Biol., 302, 205–217. algorithm using recursive variable expansion. In 2008 3rd International Rahn,R. et al. (2018) Generic accelerated sequence alignment in SeqAn using Design and Test Workshop, Monastir, 2008. IEEE, pp. vectorization and multi-threading. Bioinformatics, 34, 3437–3445. 135–140. 10.1109/IDT.2008.4802483 Rognes,T. and Seeberg,E. (2000) Six-fold speed-up of Smith–Waterman se- Landau,G.M. and Vishkin,U. (1989) Fast parallel and serial approximate quence database searches using parallel processing on common microproc- string matching. J. Algorithms, 10, 157–169. essors. Bioinformatics, 16, 699–706. Langmead,B. and Salzberg,S.L. (2012) Fast gapped-read alignment with bow- Simpson,J.T. et al. (2009) ABySS: a parallel assembler for short read sequence tie 2. Nat. Methods, 9, 357–359. data. Genome Res., 19, 1117–1123. Li,H. (2013) Aligning sequence reads, clone sequences and assembly contigs  Sosic,M. and  Sikic,M. (2017) Edlib: a C/Cþþ library for fast, exact sequence with bwa-mem. arXiv preprint arXiv: 1303.3997. Li,H. (2018) Minimap2: pairwise alignment for nucleotide sequences. alignment using edit distance. Bioinformatics, 33, 1394–1395. Bioinformatics, 34, 3094–3100. Suzuki,H. and Kasahara,M. (2017) Acceleration of nucleotide Li,I.T. et al. (2007) 160-fold acceleration of the smith-waterman algorithm semi-global alignment with adaptive banded dynamic programming. using a field programmable gate array (FPGA). BMC Bioinformatics, 8, BioRxiv, 130633. 185. Suzuki,H. and Kasahara,M. (2018) Introducing difference recurrence relations Liu,Y. et al. (2013) Cudaswþþ 3.0: accelerating smith-waterman protein for faster semi-global alignment of long sequences. BMC Bioinformatics, database search by coupling CPU and GPU SIMD instructions. BMC 19, 45. Bioinformatics, 14, 117. Szalkowski,A. et al. (2008) SWPS3 – fast multi-threaded vectorized Liu,Y. et al. (2014) SWAPHI-LS: Smith-waterman algorithm on Xeon Phi cop- Smith-Waterman for IBM cell/BE and 86/SSE2. BMC Res. Notes, 1, 107. rocessors for long DNA sequences. In 2014 IEEE International Conference Ukkonen,E. (1985) Algorithms for approximate string matching. Inf. Control, on Cluster Computing (CLUSTER). IEEE, pp. 257–265. 64, 100–118. Loving,J. et al. (2014) BitPAl: a bit-parallel, general integer-scoring sequence Wozniak,A. (1997) Using video-oriented instructions to speed up sequence alignment algorithm. Bioinformatics, 30, 3166–3173. comparison. Bioinformatics, 13, 145–150. Marco-Sola,S. et al. (2012) The gem mapper: fast, accurate and versatile align- Zhao,M. et al. (2013) SSW library: an SIMD Smith-Waterman C/Cþþ library ment by filtration. Nat. Methods, 9, 1185–1188. for use in genomic applications. PLoS One, 8, e82138.