Scrooge: A Fast and Memory-Frugal Genomic Sequence Aligner for CPUs, GPUs, and ASICs

Motivation: Pairwise sequence alignment is a very time-consuming step in common bioinformatics pipelines. Speeding up this step requires heuristics, efficient implementations and/or hardware acceleration. A promising candidate for all of the above is the recently proposed GenASM algorithm. We identify and address three inefficiencies in the GenASM algorithm: it has a high amount of data movement, a large memory footprint, and does some unnecessary work. Results: We propose Scrooge, a fast and memory-frugal genomic sequence aligner. Scrooge includes three novel algorithmic improvements which reduce the data movement, memory footprint, and the number of operations in the GenASM algorithm. We provide efficient open-source implementations of the Scrooge algorithm for CPUs and GPUs, which demonstrate the significant benefits of our algorithmic improvements. For long reads the CPU version of Scrooge achieves a 15x, 1.7x, and 1.9x speedup over KSW2, Edlib, and a CPU implementation of GenASM, respectively. The GPU version of Scrooge achieves a 4.2x 63x, 7.4x, 11x and 5.9x speedup over the CPU version of Scrooge, KSW2, Edlib, Darwin-GPU, and a GPU implementation of GenASM, respectively. We estimate an ASIC implementation of Scrooge to use 3.6x less chip area and 2.1x less power than a GenASM ASIC while maintaining the same throughput. Further, we systematically analyze the throughput and accuracy behavior of GenASM and Scrooge under a variety of configurations. As the optimal configuration of Scrooge depends on the computing platform, we make several observations that can help guide future implementations of Scrooge. Availability and implementation: https://github.com/CMU-SAFARI/Scrooge

We observe that GenASM (Senol Cali et al., 2020), a recent state-ofthe-art sequence alignment algorithm, has a large space for significantly improving its algorithm. GenASM uses only cheap bitwise operations and breaks the lower complexity bound of pairwise sequence alignment through its powerful windowing heuristic. Senol Cali et al. (2020) have already proven effectiveness of the GenASM algorithm and its accelerator implementation, thus we are motivated to further improve the GenASM algorithm and explore its potential on commodity hardware.
We identify three inefficiencies in the GenASM algorithm: (1) it has a large memory footprint due to the large size of the dynamic programming (DP) table, (2) it has a high amount of data movement between registers and memory due to frequent accesses to the DP table, and (3) it does some unnecessary work by calculating DP cells that are not useful for finding the final result. The three inefficiencies negatively impact both (1) software implementations running on commodity hardware (e.g., CPUs or GPUs) and (2) custom hardware (e.g., ASIC) implementations.
Software implementations on commodity hardware typically cannot fit all the data into fast on-chip memories (e.g., L1, scratchpad memory) due to the large memory footprint. This increases the latency and limits the bandwidth with which the DP table can be accessed. The high amount of data movement puts a high pressure on this bandwidth, which in turn limits throughput.
In contrast, custom hardware implementations can use arbitrarily large amounts of on-chip memory, but such a large on-chip memory with the high bandwidth requirement is uneconomical. For example, the hardware accelerator described in (Senol Cali et al., 2020) requires 76% and 54% of the total chip area and power consumption for the on-chip memory that stores the DP table.
The unnecessary work stems from computing cells that do not contain useful information for finding the final result. This applies to at least 25% of cells on average, as we calculate in Section 2.4.3. Doing unnecessary work affects software and hardware implementations equally, because both could use the wasted time to do useful work instead.
Our goal is to develop a fast and memory-frugal alignment algorithm by addressing the inefficiencies in the GenASM algorithm, and demonstrate its benefits with high-performance CPU and GPU implementations.
To this end we propose Scrooge 1 , which includes improvements to the GenASM algorithm based on three key ideas: These improvements (1) reduce the number of accesses to GenASM's DP table, (2) reduce the memory footprint of the DP table, and (3) eliminate unnecessary work. We experimentally demonstrate that our improvements yield significant benefits across multiple computing platforms. The CPU version of Scrooge achieves a 15×, 1.7×, and 1.9× speedup over KSW2 2 (Suzuki and Kasahara, 2018;Li, 2018), Edlib 2 (Šošić and Šikić, 2017), and a CPU implementation of GenASM, respectively. The GPU version of Scrooge achieves a 4.2×, 63×, 7.4×, 11×, and 5.9× speedup over the CPU version of Scrooge, KSW2, Edlib, Darwin-GPU 3 (Ahmed et al., 2020), and a GPU implementation of GenASM, respectively. We analytically estimate an ASIC implementation of Scrooge to use 3.6× less chip area and consume 2.1× less power compared to the prior state-of-the-art ASIC implementation of GenASM (Senol Cali et al., 2020), while maintaining the same throughput.
The contributions of this paper are as follows: • We develop three novel algorithmic improvements that are applicable to software and custom hardware implementations, collectively reducing the memory footprint by 24×, the number of memory accesses by 12×, and the number of entries of the DP • We systematically analyze the throughput and accuracy behavior of GenASM and Scrooge across a range of configurations based on real and simulated datasets for long and short reads. As the optimal configuration of Scrooge depends on the computing platform, we make several observations that can help guide future implementations of Scrooge 2 Materials and methods

Overview
The primary purpose of Scrooge is to accelerate pairwise sequence alignment through (1) a memory-frugal and efficient algorithm, and (2) optimized CPU, GPU, and ASIC implementations. Scrooge solves the approximate string matching (ASM) problem with the edit distance (Levenshtein, 1966) as the cost metric. That is, given two strings, text and pattern, Scrooge finds the minimum number of single-letter substitutions, insertions and deletions to convert text into pattern. Additionally, the sequence of edits that corresponds the edit distance is reported, which is called CIGAR string.
1 Scrooge is a name for miserly or frugal fictional characters, similar to how our proposed algorithm aims to be as resource-efficient as possible. 2 CPU-based The Scrooge algorithm is based on the GenASM algorithm (Section 2.2). Senol Cali et al. (2020) first proposed the GenASM algorithm, implemented it as an ASIC, and demonstrated GenASM's potential for very high throughput and resource efficiency. However, as we show in Section 2.3, the GenASM algorithm (1) exhibits a high bandwidth pressure, (2) exhibits a large memory footprint, and (3) does some unnecessary work. These inefficiencies limit GenASM's throughput and resource efficiency on both commodity and custom hardware. To improve the throughput and resource efficiency of GenASM implementations, it is critical to address these inefficiencies.
To this end we propose Scrooge's three novel algorithmic improvements to GenASM in Section 2.4. They are based on three key ideas: First, the DP table can be compressed by storing only the bitwise AND of multiple values (Section 2.4.1). Second, part of the DP table does not need to be stored because the traceback operation cannot reach these entries (Section 2.4.2). Third, part of the DP table can opportunistically be excluded from calculation if previous rows of the DP table already contain the information needed for finding the final result (Section 2.4.3).
In Section 3.2, we experimentally demonstrate that these improvements significantly increase performance on recent CPUs and GPUs. In Section 3.3, we explore the throughput behavior of GenASM with and without the proposed improvements across a range of configurations. We show in Section 3.4 that an ASIC implementation of Scrooge will have significantly reduced chip area and power consumption compared to the ASIC designed by Senol Cali et al. (2020), while maintaining the same throughput. In Section 3.5, we explore the accuracy behavior of GenASM and Scrooge across a range of configurations.

GenASM Algorithm
GenASM (Senol Cali et al., 2020) is based on the Bitap dynamic programming algorithm (Baeza-Yates and Gonnet, 1992;Wu and Manber, 1992) and improves Bitap's functionality, throughput and resource efficiency by (1) adding a novel traceback step to obtain the CIGAR string, (2) introducing a powerful windowing heuristic to reduce the algorithmic runtime complexity from cubic to linear, and (3) identifying an opportunity for intra-task parallelism.

Bitap Algorithm
Bitap uses only cheap bitwise operations to calculate the edit distance between two strings text and pattern. GenASM and Scrooge inherit this key advantage from Bitap. Bitap builds an (n+1)×(k+1) dynamic programming (DP) table R, where n=length(text) and k is the maximum number of edits considered. The entries of R are m-bit bitvectors, where m=length(pattern). Figure 1 shows an example of R after it is constructed by the Bitap algorithm. Theorem 1. The entries (bitvectors) of R can be interpreted as follows: In natural language, Theorem 1 states that the j-th bit of the bitvector R[i] [d] is 0 exactly if the suffix of text starting at character i and the suffix of pattern starting at character j differ by at most d edits. Following this interpretation, the first row d=d OP T that has a 0 in the first bit (j=0) of the leftmost column (i=0) indicates that the edit distance between text and pattern is d OP T . This bit is highlighted in pink in Figure 1.
Bitap (Algorithm 1) starts by pre-processing pattern into four pattern masks, one per character in the alphabet. The pattern mask for character X∈{A, C, G, T } is a bitvector of length m=length(pattern), with a 0 in the i-th bit if pattern[i]==X. See Figure 1 for an example.
Bitap populates the rightmost column (Line 5) and topmost row (Line 11) of R. The remaining entries are then calculated from their respective neighbors in the north (Line 13, insertion), north-east (Lines 14-15, deletion and substitution), and east (Line 16, match) through simple bitwise update rules. We refer to Senol Cali et al. (2020); Wu and Manber (1992) for detailed arguments on the correctness of Bitap. To follow the rest of this paper, it is sufficient to consider (1) the interpretation of R given in Theorem 1, and (2) the north-east data dependencies imposed by Algorithm 1 and shown in Figure 1.

GenASM's Traceback
For use-cases like read mapping, the pairwise sequence alignment algorithm should report the sequence of edits which achieves the edit distance, called the CIGAR string. Obtaining the CIGAR string involves tracing the optimal solution as a linear path through DP entries in their reverse construction order, this process is called traceback.  for d in 1 : k do 13: Unfortunately, the the original Bitap algorithm does not provide this capability (Baeza-Yates and Gonnet, 1992;Wu and Manber, 1992).  Figure 1) and go towards the north east of R. Whenever a 0 in one of these variables is traversed, the name of that variable is recorded as an edit (e.g., 'I' for an insertion). Second, it is sufficient to store only three out of the four variables (because S can be obtained by shifting D), saving both memory footprint and bandwidth.

GenASM's Windowing Heuristic
Bitap guarantees globally optimal edit distance, but this comes at the cost of cubic runtime complexity. GenASM reduces the runtime complexity to linear through a greedy windowing heuristic.
Bitap's runtime complexity is O(nk⌊m/word_size⌋), where word_size is the machine word size in bits, n=length(text), m=length(pattern) and k is the maximum number of edits considered. As long as the pattern length m is smaller than the machine word size, the runtime complexity is only O(nk). However, long reads can reach millions of base pairs in length (Payne et al., 2018), which are five orders of magnitude more than today's typical machine word length. To reduce the runtime complexity, GenASM proposes a greedy windowing heuristic. Instead of aligning text and pattern in one run of Bitap, the windowing heuristic runs Bitap multiple times as a subroutine in windows of size W. In each window, prefixes of size W characters of each sequence This approach has three advantages. First, instead of constructing a large table of n×m×k bits, only m W−O tables of W 3 bits must be constructed, saving memory footprint, data movement, and computation. Second, the Bitap subroutine now runs over constant-sized sequences, simplifying implementations. For example, DP entries can be statically assigned to processing elements (Senol Cali et al., 2020), and the data movement and the exact memory footprint are known at compile time. Third, the program flow (e.g., number of loop iterations per window) is entirely known at compile time, giving the compiler easy-to-digest optimization information.
The windowing strategy is greedy and heuristic, so it is possible that it could miss the optimal alignment and produce a sub-optimal one instead. The window size W and the window overlap O are tunable parameters. W can be understood as the range of solutions considered, similar to the band width (Ukkonen, 1985) in popular alignment implementations (e.g. (Šošić and Šikić, 2017;Suzuki and Kasahara, 2018;Li, 2018)). O can be understood as the coherence between windows or inverse greediness. We demonstrate in Section 3 that generally higher W and O improve accuracy, at the cost of lowering throughput, and choosing these values well is data dependent. Senol Cali et al. (2020) show that W=64 and O=24 achieve good throughput and accuracy for long and short read mapping (Senol Cali et al., 2020).

GenASM's Intra-Task Parallelism
Parallel processors, such as SIMD CPUs, GPUs or systolic array ASICs, can use multiple processing elements to either work on multiple problem instances at the same time (i.e., inter-task parallel) or cooperate on the same problem instance (i.e., intra-task parallel). Cooperating on the same problem instance is preferred because it completes more work while saving on memory footprint, as we show in Section 2.3. However, the original Bitap algorithm does not offer any such intra-task parallelism.
Senol Cali et al. (2020) identify that the DP entries within each northwest to south-east diagonal (one such diagonal is marked in red in Figure 1) do not depend on each other, hence they can be computed in parallel.

Inefficiencies in the GenASM Algorithm
We identify three inefficiencies in the GenASM algorithm: (1) it has a large amount of data movement, (2) it has a large memory footprint, and (3) it does some unnecessary work.
The combination of large amount of data movement and large memory footprint, which we quantify in Section 2.3.1 and 2.3.2 respectively, harms both software implementations running on commodity hardware, as well as custom hardware implementations. Commodity hardware (e.g., CPUs or GPUs) has a fixed amount of on-chip memory. The DP table might not fit into this on-chip memory, which introduces two inefficiencies: First, data now has to be moved a larger distance, which increases access latencies. Second, the high amount of data movement puts high pressure on memory bandwidth, which is scarce when accessing data residing off-chip. This causes the entire application to become memory bandwidth-bound, thus wasting compute resources and achieving sub-optimal performance. Custom hardware implementations (e.g. ASICs) can have as large on-chip memory as needed, but such a large and high-bandwidth on-chip memory comes at the cost of a large chip area and power consumption.
Doing unnecessary work trivially wastes runtime and energy. In Section 2.4.3, we identify the DP entries that are calculated needlessly by GenASM, and quantify how frequent they are.

Data Movement
We show that GenASM has a high amount of data movement in three steps: First, we count the number of bytes GenASM moves in a single window. Second, we calculate the operational intensity (Ofenbeck et al., 2014) from the data movement. Third, we use the roofline model (Williams et al., 2009) to visualize that GenASM's operational intensity is too low to saturate the compute resources of modern CPUs and GPUs.
We analyze a single GenASM window of size W=64 (see Section 2.2.3). Since GenASM simply repeats the calculation of a window in a loop, this analysis is representative for the entire program.
Number of bytes moved. Per window, GenASM generates W=64 intermediate bitvectors to initialize the topmost row of R (Line 11 of Algorithm 1). It also generates 4×64×64 intermediate bitvectors for the inner entries of R (Lines 13-17 of Algorithm 1), out of which 3×64×64 are stored for traceback (Section 2.2.2). Bitvectors have size W=64 bits each. Assuming the processing elements communicate neighbor entry values without memory accesses 4 , each of these bitvectors is written to memory exactly once per window (during the construction of R), and not read again 5 . This yields data_movement = (64 + 3×64×64) × 64b = 98,816B = 96.5KiB.
Operational intensity. The operational intensity is defined as the ratio of work done over data movement (Ofenbeck et al., 2014). We define work as the number of arithmetic and logic 64-bit integer operations (op), and data_movement as the number of bytes moved between registers and the memory hierarchy and scratchpad memory.
The GenASM algorithm's operational intensity is then I = work data_movement ≈ 0.29 op B . Roofline model. The roofline model (Williams et al., 2009) plots the upper limit of achievable compute throughput for different operational intensities for a given processor. It consists of horizontal peak compute throughput rooflines, and sloped memory bandwidth rooflines. The intersection point between the operational intensity and the roofline (either compute or memory) gives an upper bound on the achievable compute throughput. We use the roofline model to argue that GenASM's implementation on CPUs and GPUs exhibits a high amount of data movement relative to the respectively available bandwidths. Figure 2 shows the roofline plots for an NVIDIA A6000 GPU 6 (NVIDIA, 2020) and an Intel Xeon Gold 5118 CPU 7 (Intel, 2017), including their respective on-chip scratchpad 8 , cache, and offchip 9 memory bandwidths (drawn in shades of blue) and peak compute throughputs (draw in shades of green). GenASM's operational intensity is drawn in red.
From Figure 2 we make three observations. First, GenASM is heavily memory bandwidth-bound for a modern CPU and GPU if the data resides off-chip (the red and dark blue lines intersect far below the green line). Second, GenASM would no longer be memory bandwidth-bound if it moves data off-chip over an order of magnitude less frequently (the red line has to be shifted far to the right for the red-dark blue intersection point to reach the green line). Third, if the data resides in a fast on-chip memory instead, GenASM can reach peak compute throughput, even with the high amount of data movement in the baseline algorithm (the red and light blue lines intersect above the green line). However, as we show in Section 2.3.2, the memory footprint is too large to fit into the available onchip memory for some commodity hardware, and building large enough on-chip memories it uneconomical.
Based on these observations, we conclude that (1) GenASM cannot saturate commodity hardware with computation, and (2) data movement should be reduced to address this inefficiency.

Memory Footprint
In this section, we demonstrate the overheads associated with GenASM's large memory footprint.
GenASM keeps all written bitvectors in memory for the traceback step. Thus, directly following from the calculated data movement in Section 2.3.1, the memory footprint is 96.5KiB. For comparison, NVIDIA's most recent GPU microarchitecture, Ampere, provides up to 99KiB of high-bandwidth on-chip scratchpad memory per GPU core (streaming multiprocessor, SM in NVIDIA terminology) (NVIDIA, 2021b). Thus, one SM can hold the DP table for exactly one GenASM problem instance in its scratchpad. One thread block of two warps 6 We obtain the global memory bandwidth from the datasheet. We derive compute throughput and shared memory bandwidth from the CUDA Programming Guide (NVIDIA, 2021b) by multiplying the respective per-cycle throughputs with the frequency and number of streaming multiprocessors (SMs). 7 We obtain the cache bandwidths and compute throughputs from the Intel Optimization Manual (Intel, 2022) by multiplying the respective per-cycle throughputs with frequency and number of cores. We obtain the off-chip (DRAM) bandwidth by assuming 2,400MHz DDR4 memory, multiplied by the maximum number of memory channels supported by the CPU. 8 called shared memory in CUDA 9 called global memory in CUDA, DRAM for CPUs (i.e., 2×32 threads) can work on single GenASM problem instance, however this does not saturate the compute resources in the SM. This is because modern GPUs are designed to alternate between executing multiple independent instruction streams for the purpose of hiding the latency of instructions (Arafa et al., 2019). The under-saturation of the compute resources in the SM due to too few independent instruction streams is called low occupancy and causes the unused computational resources to be wasted (NVIDIA, 2021b). Hence, the occupancy should be increased by working on multiple problem instances per SM. Multiple problem instances can fit into memory by either reducing the memory footprint per problem instance, or placing the DP tables into the GPU's off-chip memory, which has a much larger capacity. However, we show in Section 2.3.1 that the latter is not an efficient a solution due to the off-chip memory's limited bandwidth.
Custom hardware implementations (e.g. ASICs) can have as large onchip memory as needed. For example, the GenASM ASIC designed by Senol Cali et al. (2020) uses scratchpads of 96.5 KiB each to hold the DP tables. However, these scratchpads occupy 76% of the total chip area and consume over 54% of the chip power. This limits the performance achievable with a given chip area and power budget.
In summary, GenASM has a large memory footprint compared to typical on-chip memory capacities in commodity hardware, and while sufficiently large on-chip scratchpads can be designed for custom hardware implementations, it is uneconomical to do so.

Scrooge
We have shown in Section 2.3 that GenASM has a high amount of data movement and high memory footprint per problem instance. We have elaborated that this combination either limits performance (on commodity hardware), or requires expensive large caches (on custom hardware), both of which are undesirable. Thus, our strategy is to reduce the GenASM algorithm's memory footprint as much as possible while introducing minimal computational overhead. We present three novel algorithmic improvements that collectively achieve a 24× reduction in memory footprint, as well as a 12× reduction in data movement.

Improvement 1 -Store Entries, not Edges (SENE)
As we explain in Section 2.2.2, GenASM stores 3 bitvectors per entry of the table R to enable traceback. If we imagine a graph where the entries of R are nodes and the intermediate bitvectors are edges connecting their source and target entries, GenASM stores 3 ingoing edges for most nodes (Figure 3). We propose SENE, observing that it is possible to regenerate (only) the required edges on-demand during traceback, if the nodes (entries of table R) are stored instead. Regenerating edges is done by applying the update rules in Algorithm 1 on requested neighbors entries during traceback.  Due to the fixed word sizes and word alignment requirements of commodity hardware, the number of bits stored for each bitvector cannot be chosen freely, and we observe GenASM's choice of O=24 does not yield the best performance. For example, we show in Section 3.2 that for a modern GPU O = 33 achieves the best throughput results for W = 64, because the stored bitvectors perfectly fit into a 32-bit word. Note that increasing O improves accuracy, see Section 2.2.3 for an intuition and Section 3.5 for measurements.
Cost and benefits. First, determining which bitvectors to store and extracting the left half of each stored bitvector incurs a small computational overhead. Second, increasing O from 24 to 33 means the algorithm makes 9 characters less progress per window.
By discarding the right half of each bitvector and the rightmost O columns of R, DENT improves the memory footprint by We describe in the supplementary materials how DENT can be extended to rows for a total 8× memory footprint reduction.
By the same calculation, the number of writes to table R is reduced by approximately 4×, assuming the forefront diagonal (marked red in Figure 1) is kept in registers and communicated directly 4 . Fig. 4. DENT exploits that the windowing heuristic stops traceback before the entire path is traversed. The area never reached by traceback can be discarded.

Improvement 3 -Early Termination (ET)
The edit distance is determined by the highest row of R that contains a 0 in the most significant bit in the leftmost column. Traceback is also started from this entry, and any traceback path can only go to the right and up. We observe that at no point do the rows of higher cost than distance(pattern, text) contain useful information (See Figure 5 for an example). Thus, we propose building R row-wise, and terminating the algorithm early as soon as the most significant bit in the first entry of the current row is a 0. Cost and benefits. Early Termination does not yield a constant factor improvement in either memory footprint or runtime: If distance(pattern, text) = W we are not able to terminate early at all. However, typical input pairs incur fewer edits. For correct candidate pairs, the edit distance will be low, e.g. up to 15% for long reads (Alser et al., 2021). Even uncorrelated random sequence pairs of length W over a 4letter alphabet have a Hamming distance (Hamming, 1950) of only 3 4 W on average, and the edit distance is at most the Hamming distance. Thus, on average Scrooge can skip at least 25% of the entries of R, saving computation as well as data movement. The early termination improvement is required for the extended DENT improvement, as we describe in the supplementary materials.
Conflict with intra-task parallelism. Recall that GenASM provides the option for intra-task parallelism, as presented in Section 2.2.4. Exploiting this parallelism requires the available processing elements to build R in a diagonal-wise fashion, as shown in Figure 1. However, as we describe in Section 2.4.3, to make full use of Early Termination, R should be built row-wise. As a compromise, we implement Early Termination in diagonal-wise fashion in our GPU implementation. As in the row-wise version, construction on R stops as soon as the leftmost processing element finds a 0 in the most significant bit. Due to the diagonal-wise computation, the other processing elements have already computed several rows ahead at this point, i.e. done unnecessary work. For this reason the benefit of Early Termination is limited in intra-task parallel implementations, such as our GPU implementation, while being much more significant in rowwise implementations, such as our CPU implementation. We reaffirm these effects experimentally in Section 3.3.

Implementation
We implement C++ versions of our algorithm for CPUs and NVIDIA GPUs. They are exposed as simple library functions for pairwise alignment. Each improvement and implementation constant can be easily configured at compile time through preprocessor macros. The implementations, as well as baselines and evaluation scripts are available at https://github.com/cmu-safari/Scrooge.
CPU. The CPU version converts the input pairs to a two-bit-perbasepair encoding, but padded to 8 bits. Each thread works on a single pairwise alignment at a time and obtains sequence pairs from a global queue with workstealing. During each call to the Bitap subroutine, the thread calculates the DP table R in a row-wise fashion.
GPU. The GPU version is implemented using CUDA 11.1 (NVIDIA, 2021b) and targets GPUs of compute capability 7.0 and higher (NVIDIA, 2021a). The input sequence pairs are converted to a two-bit-per-basepair encoding and transferred to the GPU. Each threadblock works on a single pairwise alignment at a time and obtains sequence pairs from a global queue with workstealing. During each call to the Bitap subroutine, the threadblock calculates the DP table R in a diagonal-wise fashion, and each of the W threads in the threadblock calculates a single column of R. Threads resolve their mutual data dependencies using warp shuffle instructions within a warp, and using shared memory across warps. A single thread per warp executes the traceback operation. The size of the CIGAR string is not known ahead of time, hence it is stored as a linked list in global memory.
We evaluate the throughput and accuracy of Scrooge in three classes of experiments: First, we compare the throughputs of all evaluated tools and show that Scrooge outperforms state-of-the-art aligners. Second, we evaluate the throughput benefits of Scrooge's algorithmic improvements, and its sensitivity to different choices for W an O. Third, we evaluate Scrooge's accuracy compared to Edlib. We define throughput as the number of pairwise alignments per second for a given dataset and accuracy as the similarity between the alignment scores of Scrooge's output and the output of Edlib (Šošić and Šikić, 2017).
We run all CPU evaluations on a dual socket Intel Xeon Gold 5118 (2× 12 physical cores, 2× 24 logical cores) at 3.2GHz with 196GiB DDR4 RAM. We run all GPU evaluations on an NVIDIA A6000 (NVIDIA, 2020).

Datasets
We simulate 115,240 PacBio reads from the human genome using PBSIM2 (Ono et al., 2020), each of length 10 kilo bases and with a target error rate of 5%. We obtain the ground truth locations (one per read) from PBSIM2, thus obtaining 115,240 candidate pairs for our PBSIM2 groundtruth dataset. We map 500 of the simlated PacBio reads to the human genome using minimap2 (Li, 2018) and obtain all chains (candidate locations) it generates using the -P flag, 138,929 locations in total. This constitutes our PBSIM2 chained dataset. We map 100,000 Illumina short reads from the dataset with accession number SRR13278681 to the human genome using minimap2 (Li, 2018) and obtain all chains (candidate locations) it generates using the -P flag, 9,612,222 locations in total. This constitutes our Illumina chained dataset. More statistics about the datasets can be found in the supplementary materials. The exact datasets and commandlines that produced them are available in our github repository.

Throughput
We run the CPU-based tools using 48 threads, and take the fastest configuration from a parameter sweep for Darwin-GPU 10 . We empirically configure Scrooge's CPU and GPU implementations with W=64, O=33 for the PBSIM2 chained dataset and W=32, O=17 for the Illumina chained dataset, and enable the combinations of improvements that yield the best throughput. Figure 6 shows Scrooge's algorithmic improvements significantly speed up the alignment of long and short reads over all baselines. In particular, the CPU implementation of Scrooge has a 1.9× higher throughput than our CPU implementation of GenASM for long reads and and a 4.2× higher throughput for short reads. The GPU implementation of Scrooge has a 5.9× higher throughput than our GPU implementation of GenASM for long reads and a 2.3× higher throughput for short reads. The CPU and GPU speedups over GenASM are entirely due to Scrooge's algorithmic improvements, since our Scrooge and GenASM implementations are similarly optimized.
Note that WFA, KSW2, CUDASW++3.0 and Darwin solve a more general formulation of the alignment problem with affine gap scores (Gotoh, 1982). This puts them at a performance disadvantage. In contrast, Edlib (Šošić and Šikić, 2017), GenASM (Senol Cali et al., 2020), and Scrooge solve a less general but more efficient formulation of the alignment problem with unit costs (edit distance or Levenshtein distance (Levenshtein, 1966)).  Fig. 6. Scrooge's throughput relative to CPU and GPU baselines for the affine gap cost (KSW2, Darwin-GPU) and edit distance metrics (Edlib) for long and short reads.

Sensitivity Analysis
We explore the throughput benefits of our algorithmic improvements in parameter sweeps over (1) the number of GPU and CPU threads, (2) the window size (W) parameter, and (3) the window overlap (O) parameter. GPU threads. First, we run a scaling experiment on a GPU for GenASM, Scrooge with the SENE improvement, Scrooge with the DENT improvement, and Scrooge with all three proposed improvements. Based on Figure 7, we make five observations: First, we observe that SENE and DENT individually improve performance when the DP table is placed in either shared or global memory. Second, we observe that SENE, DENT and Early Termination can be combined for greater benefits. Third, we observe that shared memory achieves the best performance, but only when all proposed memory improvements are applied. Fourth, in some configurations, using global (off-chip) memory can be faster and shared (on-chip) memory. This is because of the low occupancy due to the small shared memory capacity (Section 2.3.2), while global memory can fit many problem instances and thus allows for higher occupancy. Throughput then scales until eventually the memory bandwidth to global memory limits performance. Finally, we observe that the baseline GenASM algorithm cannot run using shared memory at all, although we showed in Section 2.3.2 that a single instance of the baseline DP table has a footprint of 98.5KiB and thus should use fit into the 100KiB of shared memory. This is because our implementation requires some additional memory, such as for communication between processing elements. Thus, we cannot fit even a single instance into shared memory.
We also evaluate the performance benefits of each of the 2 3 = 8 possible combinations of the SENE, DENT and Early Termination improvements. From Supplementary Figure 8 we observe that the Early Termination improvement consistently improves performance of the GPU implementation across all combinations, but only by a few percent.
CPU threads. We run a similar scaling experiment on a CPU for GenASM, Scrooge with the SENE improvement, Scrooge with the Early Termination (ET) improvement, and Scrooge with the SENE and Early Termination improvements. From Figure 8a we make three observations: First, we observe that Early Termination improves 10 For a meaningful comparison we ensure that Darwin-GPU's alignment component fully aligns all sequence pairs. We disable conflicting heuristics in configuration files where possible, and modify the code where necessary. performance significantly. This is in contrast to our GPU implementation, where Early Termination did not show significant benefits. This is because our CPU implementation builds the DP table R row-wise, while our GPU implementation builds R diagonal-wise (see Section 2.4.3). Second, SENE improves performance consistently, but less significantly than in the GPU case. This is because modern CPUs have relatively large onchip cache capacities (e.g., 1MiB L2 cache per core on the Xeon Gold 5118 we evaluated on (Intel, 2017)). Thus, the DP table easily fits into the L2 cache even without Scrooge's algorithmic improvements, and hence reducing the memory footprint is not as important. Third, Scrooge scales linearly up to 24 threads, but does not scale at all from 24 to 48 threads. The system we evaluate on has 24 physical cores and 48 logical cores, thus we hypothesize that our CPU implementation does not benefit from simultaneous multithreading (Hyper-Threading in Intel terminology) (Marr et al., 2002) due to the low latencies of bitwise instructions (Fog, 2021), which is what the majority of our implementation consists of. We also evaluate the performance benefits of each of the 2 3 =8 possible combinations of the SENE, DENT and Early Termination improvements. From Supplementary Figure 9a we observe that DENT does not yield performance benefits on a CPU, but actually reduces performance. This is because of two reasons: First, the memory footprint reduction is not critical on modern CPUs, as we explained for SENE. Second, although extracting the right bits for DENT to store (Section 2.4.2) is conceptually simple, it requires several additional instructions to be executed per bitvector.
The three key takeaways from these experiments are that (1) the SENE and DENT memory improvements yield significant benefits if performance is limited by memory bandwidth or capacity, (2) some of the algorithmic improvements can cause a performance loss in practice, and (3) the ideal combination of improvements depends on the computation platform (e.g., the available on-chip cache capacity) and the exact implementation (e.g., row-wise or diagonal-wise). Window size (W). We explore the sensitivity of Scrooge's throughput to the window size parameter W (Section 2.2.3) on CPUs and GPUs. We vary W and set O=W//2+1. Note that larger W improve accuracy (Section 2.2.3).
From the CPU results in Figure 8b we make two observations: First, we observe that performance generally reduces as W increases. This is because the number of calculated bits per window increases cubically with increasing W. Second, we observe a sudden throughput dropoff when W increases past 64. This is because the word size of the Xeon Gold 5118 CPU is 64 bits, thus if W>64 each bitvector operation has to be emulated using multiple word-sized machine instructions. This emulation is conceptually simple (e.g. carry over shifted bits), but requires several additional instructions, causing the performance dropoff. We show the results for each of the 2 3 = 8 possible combinations of the SENE, DENT and Early Termination improvements in Supplementary Figure 9b.
From the GPU results in Figure 9 we make five observations: First, as on the CPU, performance generally reduces as W increases. Second, as on the CPU, there are sudden performance dropoffs whenever the window (and thus bitvector) size surpasses a multiple of the machine word size (32 bits on the tested GPU (NVIDIA, 2021b)). Third, some configurations can run only for small W when the DP table R is placed in shared memory. This is because as W increases, the memory footprint of R increases cubically, and some configurations run out of shared memory. With global memory, the DP tables can occupy as much of the GPUs off-chip memory as needed, thus Scrooge can run even very large W, although the performance eventually becomes impractically low. Fourth, for very small W, shared memory can offer significantly better performance than global memory. This is because for smaller W the memory footprint per DP table is smaller, and more instances can be kept in shared memory, thus easily achieving high occupancy. In contrast, global memory remains bound by memory bandwidth, because the operational intensity is unchanged. Fifth, for small window sizes and shared memory it is beneficial to run without the memory improvements SENE an DENT. This corroborates the observation made for the CPU implementation: If memory capacity and bandwidth are not limiting factors, then the computational overheads of SENE and DENT overshadow their benefits. In contrast, when W is large, the memory improvements become increasingly important because of the inverse argument. We show the results for each of the 2 3 = 8 possible combinations of the SENE, DENT and Early Termination improvements in Supplementary Figure 10.
The three key takeaways from these experiments are that (1) the ideal combination of improvements depends on the window size (W) parameter, (2) the memory improvements SENE and DENT memory improvements are increasingly beneficial as the window size increases, and (3)  From the GPU results in Figure 9 we make two observations: First, when DENT is disabled, performance decreases strictly and smoothly as O increases. This is because the algorithm makes W-O characters progress per window, i.e. progress per window decreases linearly with increasing O. Second, we observe that with DENT enabled, performance can sometimes increase as O increases, up to some limit. This is because DENT reduces the memory footprint of the DP table more as O increases. We observe a large and sudden increase in performance at O=34. This is because with DENT, each bitvector stored for traceback has size W-O+1, thus for O≥33 each stored bitvector fits into a single 32-bit machine word (NVIDIA, 2020). We show the results for each of the 2 3 = 8 possible combinations of the SENE, DENT and Early Termination improvements in Supplementary  Figure 11.
The key takeaway from this experiment is that a larger window overlap (O) can improve performance up to some limit. This is convenient because larger O also improve accuracy. In particular, we found for W=64 the optimal performance is attained by O=33, increasing both throughput and accuracy over the default operating point chosen by Senol Cali et al.

Benefits for hardware implementations
The ASIC designed by Senol Cali et al. (2020) uses a large on-chip scratchpad to store bitvectors for traceback. This scratchpad alone accounts for 0.256mm 2 (76%) of silicon area and 0.055W (54%) of power out of a total of 0.334mm 2 and 0.101W per accelerator core. Our proposed algorithmic improvements can be applied to their design through minor modifications. We estimate the area and power cost of such an ASIC implementation of Scrooge analytically as follows: 1. We start with the DC-logic, DC-SRAM and TB-logic area and power numbers reported in (Senol Cali et al., 2020) 2. We estimate Scrooge's TB-SRAM area and power cost with CACTI 7 (Balasubramonian et al., 2017), as in (Senol Cali et al., 2020), but with Scrooge's reduced memory footprint and data movement numbers. 3. We account for the logic overhead of SENE by adding the area and power of a single DC processing element (Senol Cali et al., 2020) to the traceback (TB) logic cost, this accounts for recomputing edges during traceback. We assume no overhead for SENE during the construction of R, since the ANDed bitvectors are already computed. 4. We assume no overheads for applying DENT, since it simply masks out bits when storing the bitvectors, which is trivial in hardware. Table 1 lists the area and power breakdowns obtained with this methodology, and the breakdown of (Senol Cali et al., 2020) as a comparison point. In particular, we observe a 3.6× reduction in chip area and a 2.1× reduction in chip power consumption, while maintaining the same throughput. These improvements come from (1) the reduced TB SRAM capacity, and (2) the reduced TB SRAM bandwidth.
The key takeaway from this estimate is that Scrooge's algorithmic improvements (1) are directly applicable to and (2) yield significant benefits for an ASIC implementation of GenASM.

Accuracy
The GenASM algorithm (Senol Cali et al., 2020), which Scrooge is based on, is a greedy heuristic algorithm, as explained in Section 2.2. Our improvements do not introduce additional inaccuracy, in fact Scrooge's default operating point of W=64 O=33 increases accuracy over GenASM's default operating point of W=64 O=24 (Senol Cali et al., 2020). Thus the following analysis explores the accuracy of the GenASM algorithm itself at different operating points. We analyze the accuracy compared to optimal edit distance solutions, such as Edlib (Šošić and Šikić, 2017). We evaluate the generated alignments based on minimap2's default affine gap scoring model. We run two types of experiments: First, we sweep over the window size W. Second, we sweep over the window overlap W. Window size (W). We explore the sensitivity of Scrooge's accuracy to the window size parameter W (Section 2.2.3). We vary W and set O=W//2+1. For each experiment, we record the 0.5, 0.1, 0.01 and 0.001 percentile alignment scores 11 of Scrooge/GenASM and compare to Edlib as an ideal upper bound.
From Figure 11 we make three observations: First, the accuracy depends on the dataset. Second, small window sizes are sufficient for GenASM to find the optimal edit distance alignment for most of the sequence pairs. For example, the median alignment score is already optimal at W=32 for the PBSIM2 groundtruth dataset, and at W=8 for the Illumina chained dataset. Third, to find the optimal alignment for a few worst-case pairs large window sizes are required: For example, the optimal alignment for the 0.001 percentile in the PBSIM2 groundtruth dataset is only found for W≥80. We manually inspect several of these "difficult" sequence pairs to find the reason for their apparent difficulty. We observe sequence pairs are aligned poorly when they contain extremely noisy and repetitive sub-sequences. However, these pairs will be aligned optimally if the window size is larger than the noisy sub-sequences. We provide a tool for inspecting such noisy sequence pairs and the alignment they generate. Window overlap (O). We explore the sensitivity of Scrooge's accuracy to the window overlap parameter O (Section 2.2.3). We sweep over O and run experiments for each W ∈ {32, 64, 96, 128}. For each experiment we record the 0.01 percentile alignment score of Scrooge/GenASM and compare to Edlib as an ideal upper bound.
From Figure 12 we make two observations: First, accuracy improves as O increases. Second, we observe that W and O need to be balanced to achieve good accuracy. For example, choosing too small W=32 for the PBSIM2 groundtruth dataset cannot be made up with even large O=30. Similarly, choosing O close to 0 hurts accuracy for both datasets, even when W is large.
The two key takeaways from these experiments are that (1) W and O need to be chosen per dataset, and (2) W and O should be balanced, i.e., increased or reduced together for the best accuracy.

Discussion
To our knowledge this is the first paper to explore the benefits of the GenASM algorithm on commodity hardware and its accuracy under different configurations. We have systematically analyzed the inefficiencies in the GenASM algorithm and proposed three novel improvements to address them, and have demonstrated the corresponding benefits experimentally for commodity hardware, and analytically for custom hardware.
We have already extensively compared to WFA (Marco-Sola et al., 2020), KSW2 (Suzuki and Kasahara, 2018;Li, 2018), Edlib (Šošić and Šikić, 2017), CUDASW++3.0 (Liu et al., 2013), and Darwin-GPU (Ahmed et al., 2020). There are several other works that accelerate sequence alignment: NVBIO (NVIDIA, 2014) is a multi-purpose library for accelerating bioinformatics applications using GPUs, but is no longer maintained. Gasal2 (Ahmed et al., 2019) is a recent GPU aligner limited to short reads. CUDAlign4.0 (de Oliveira Sandes et al., 2016) can efficiently align a single pair of extremely long (chromosome sized) sequences, with use cases such as whole genome alignment. Adept (Awan et al., 2020) is a recent GPU aligner for short and long reads but does not support traceback, i.e. only reports the alignment score.
The Darwin accelerator (Turakhia et al., 2019) implements a Smith-Waterman-Gotoh accelerator for long reads using a similar greedy strategy to GenASM called tiling. We have compared to the GPU implementation of this algorithm, Darwin-GPU. With the feasibility of such greedy algorithms already demonstrated, a possible future direction is to study the suitability of different algorithms to such greedy heuristics, such as Myers' bitvector algorithm (Myers, 1999) or the recently proposed wavefront algorithm (Marco-Sola et al., 2020). Another interesting direction will be to explore the effectiveness of our improvements for other implementations of greedy windowing or tiling, like Darwin. We believe the DENT improvement can be applied directly to Darwin.