Genome sequencing, which determines the DNA sequence of an organism, plays a pivotal role in enabling many medical and scienti c advancements in personalized medicine [6, 20, 34, 53, 59] , evolutionary theory [46, 139, 140] , and forensics [17, 25, 179] . Modern genome sequencing machines [77-79, 132-135, 152] can rapidly generate massive amounts of genomics data at low cost [8, 118, 153] , but are unable to extract an organism's complete DNA in one piece. Instead, these machines extract smaller random fragments of the original DNA sequence, known as reads. These reads then pass through a computational process known as read mapping, which takes each read, aligns it to one or more possible locations within the reference genome, and nds the matches and di erences (i.e., distance) between the read and the reference genome segment at that location [6, 177] . Read mapping is the rst key step in genome sequence analysis.
State-of-the-art sequencing machines produce broadly one of two kinds of reads. Short reads (consisting of no more than a few hundred DNA base pairs [30, 158] ) are generated using short-read sequencing (SRS) technologies [144, 164] , which have been on the market for more than a decade. Because each read fragment is so short compared to the entire DNA (e.g., a human's DNA consists of over 3 billion base pairs [166] ), short reads incur a number of reproducibility (e.g., non-deterministic mapping) and computational challenges [7, 10, 12, 52, 118, 159, [176] [177] [178] . Long reads (consisting of thousands to millions of DNA base pairs) are generated using long-read sequencing (LRS) technologies, of which Oxford Nanopore Technologies' (ONT) nanopore sequencing [26, 35, 40, 82, 83, 89, 97, 112, 113, 116, 143, 152] and Paci c Biosciences' (PacBio) single-molecule real-time (SMRT) sequencing [18, 47, 114, 123, 145, 146, 165, 171] are the most widely used ones. LRS technologies are relatively new, and they avoid many of the challenges faced by short reads.
LRS technologies have three key advantages compared to SRS technologies. First, LRS devices can generate very long reads, which (1) reduces the non-deterministic mapping problem faced by short reads, as long reads are signi cantly more likely to be unique and therefore have fewer potential mapping locations in the reference genome; and (2) span larger parts of the repeated or complex regions of a genome, enabling detection of genetic variations that might exist in these regions [165] . Second, LRS devices perform real-time sequencing, and can enable concurrent sequencing and analysis [111, 142, 146] . Third, ONT's pocket-sized device (Min-ION [133] ) provides portability, making sequencing possible at remote places using laptops or mobile devices. This enables a number of new applications, such as rapid infection diagnosis and outbreak tracing (e.g., COVID-19, Ebola, Zika, swine u [37, 48, 64, 68, 85, 142, 167, 173] ). Unfortunately, LRS devices are much more error-prone in sequencing (with a typical error rate of 10-15% [19, 83, 165, 170] ) compared to SRS devices (typically 0.1% [60, 61, 141] ), which leads to new computational challenges [152] .
For both short and long reads, multiple steps of read mapping must account for the sequencing errors, and for the differences caused by genetic mutations and variations. These errors and di erences take the form of base insertions, deletions, and/or substitutions [121, 125, 154, 163, 169, 174] . As a result, read mapping must perform approximate (or fuzzy) string matching (ASM). Several algorithms exist for ASM, but stateof-the-art read mapping tools typically make use of an expen-sive dynamic programming based algorithm [100, 126, 154] that scales quadratically in both execution time and required storage. This ASM algorithm has been shown to be the major bottleneck in read mapping [8, 10, 55, 66, 75, 122, 162] . Unfortunately, as sequencing technologies advance, the growth in the rate that sequencing devices generate reads is far outpacing the corresponding growth in computational power [8, 32] , placing greater pressure on the ASM bottleneck. Beyond read mapping, ASM is a key technique for other bioinformatics problems such as whole genome alignment (WGA) [27, 28, 41, 42, 70, 95, 102, 106, 115, 151, 160] and multiple sequence alignment (MSA) [29, 45, 69, 98, 107, 127, 128, 136, 150] , where two or more whole genomes, or regions of multiple genomes (from the same or di erent species), are compared to determine their similarity for predicting evolutionary relationships or nding common regions (e.g., genes). Thus, there is a pressing need to develop techniques for genome sequence analysis that provide fast and e cient ASM.
In this work, we propose GenASM, an ASM acceleration framework for genome sequence analysis. Our goal is to design a fast, e cient, and exible framework for both short and long reads, which can be used to accelerate multiple steps of the genome sequence analysis pipeline. To avoid implementing more complex hardware for the dynamic programming based algorithm [22, 33, 49, 65, 87, 88, 147, 162] , we base GenASM upon the Bitap algorithm [21, 174] . Bitap uses only fast and simple bitwise operations to perform approximate string matching, making it amenable to e cient hardware acceleration. To our knowledge, GenASM is the rst work that enhances and accelerates Bitap.
To use Bitap for GenASM, we make two key algorithmic modi cations that allow us to overcome key limitations that prevent the original Bitap algorithm from being e cient for genome sequence analysis (we discuss these limitations in Section 2.3). First, to improve Bitap's applicability to di erent sequencing technologies and its performance, we (1) modify the algorithm to support long reads (in addition to already supported short reads), and (2) eliminate loop-carried data dependencies so that we can parallelize a single string matching operation. Second, we develop a novel Bitap-compatible algorithm for traceback, a method that utilizes information collected during ASM about the di erent types of errors to identify the optimal alignment of reads. The original Bitap algorithm is not capable of performing traceback.
In GenASM, we co-design our modi ed Bitap algorithm and our new Bitap-compatible traceback algorithm with an areaand power-e cient hardware accelerator, which consists of two components: (1) GenASM-DC, which provides hardware support to e ciently execute our modi ed Bitap algorithm to generate bitvectors (each of which represents one of the four possible cases: match, insertion, deletion, or substitution) and perform distance calculation (DC) (which calculates the minimum number of errors between the read and the reference segment); and (2) GenASM-TB, which provides hardware support to e ciently execute our novel traceback (TB) algorithm to nd the optimal alignment of a read, using the bitvectors generated by GenASM-DC. Our hardware accelerator (1) balances the compute resources with available memory capacity and bandwidth per compute unit to avoid wasting resources, (2) achieves high performance and power e ciency by using specialized compute units that we design to exploit data locality, and (3) scales linearly in performance with the number of parallel compute units that we add to the system. Use Cases. GenASM is an e cient framework for accelerating genome sequence analysis that has multiple possible use cases. In this paper, we describe and rigorously evaluate three use cases of GenASM. First, we show that GenASM can e ectively accelerate the read alignment step of read mapping (Section 10.2). Second, we illustrate that GenASM can be employed as the most e cient (to date) pre-alignment lter [9, 10] for short reads (Section 10.3). Third, we demonstrate how GenASM can e ciently nd the edit distance (i.e., Levenshtein distance [100] ) between two sequences of arbitrary lengths (Section 10.4). In addition, GenASM can be utilized in several other parts of genome sequence analysis as well as in text analysis, which we brie y discuss in Section 11.
Results Summary. We evaluate GenASM for three different use cases of ASM in genome sequence analysis using a combination of the synthesized SystemVerilog model of our hardware accelerators and detailed simulation-based performance modeling. (1) For read alignment, we compare GenASM to state-of-the-art software (Minimap2 [102] and BWA-MEM [101] ) and hardware approaches (GACT in Darwin [162] and SillaX in GenAx [55] ), and nd that GenASM is signi cantly more e cient in terms of both speed and power consumption. For this use case, we compare GenASM only with the read alignment steps of the baseline tools and accelerators. For long reads, GenASM achieves 116× and 648× speedup over 12-thread runs of the alignment steps of Min-imap2 and BWA-MEM, respectively, while reducing power consumption by 37× and 34×. Compared to GACT, GenASM provides 6.6× the throughput per unit area and 10.5× the throughput per unit power for long reads. For short reads, GenASM achieves 158× and 111× speedup over 12-thread runs of the alignment steps of Minimap2 and BWA-MEM, respectively, while reducing power consumption by 31× and 33×. Compared to SillaX, GenASM is 1.9× faster at a comparable area and power consumption. (2) For pre-alignment ltering of short reads, we compare GenASM with a state-ofthe-art FPGA-based lter, Shouji [9] . GenASM provides 3.7× speedup over Shouji, while reducing power consumption by 1.7×, and also signi cantly improving the ltering accuracy.
(3) For edit distance calculation, we compare GenASM with a state-of-the-art software library, Edlib [155] , and FPGAbased accelerator, ASAP [22] . Compared to Edlib, GenASM provides 22-12501× speedup, for varying sequence lengths and similarity values, while reducing power consumption by 548-582×. Compared to ASAP, GenASM provides 9.3-400× speedup, while reducing power consumption by 67×. This paper makes the following contributions: • To our knowledge, GenASM is the rst work that enhances and accelerates the Bitap algorithm for approximate string matching. We modify Bitap to add e cient support for long reads and enable parallelism within each ASM operation. We also propose the rst Bitap-compatible traceback algorithm. We open source our software implementations of the GenASM algorithms [148] . • We present GenASM, a novel approximate string matching acceleration framework for genome sequence analysis. GenASM is a power-and area-e cient hardware implementation of our new Bitap-based algorithms. • We show that GenASM can accelerate three use cases of approximate string matching (ASM) in genome sequence analysis (i.e., read alignment, pre-alignment ltering, edit distance calculation). We nd that GenASM is greatly faster and more power-e cient for all three use cases than stateof-the-art software and hardware baselines.
A common approach to the rst step in genome sequence analysis is to perform read mapping, where each read of an organism's sequenced genome is matched against the reference genome for the organism's species to nd the read's original location. As Figure 1 shows, typical read mapping [6, 96, 101, 102, 105, 177] is a four-step process. First, read mapping starts with indexing 0 , which is an o ine preprocessing step performed on a known reference genome. Second, once a sequencing machine generates reads from a DNA sequence, the seeding process 1 queries the index structure to determine the candidate (i.e., potential) mapping locations of each read in the reference genome using substrings (i.e., seeds) from each read. Third, for each read, pre-alignment ltering 2 uses ltering heuristics to examine the similarity between a read and the portion of the reference genome at each of the read's candidate mapping locations. These ltering heuristics aim to eliminate most of the dissimilar pairs of reads and candidate mapping locations to decrease the number of required alignments in the next step. Fourth, for all of the remaining candidate mapping locations, read alignment 3 runs a dynamic programming based algorithm to determine which of the candidate mapping locations in the reference matches best with the input read. As part of this step, traceback is performed between the reference and the input read to nd the optimal alignment, which is the alignment with the highest likelihood of being correct (based on a scoring function [62, 117, 168] ). The optimal alignment is de ned using a CIGAR string [103] , which shows the sequence and position of each match, substitution, insertion, and deletion for the read with respect to the selected mapping location of the reference.
Hash- 
The goal of approximate string matching [125] is to detect the di erences and similarities between two sequences. Given a query read sequence Q=[q 1 q 2 . . . q m ], a reference text sequence T =[t 1 t 2 . . . t n ] (where m = |Q|, n = |T |, n ≥ m), and an edit distance threshold E, the approximate string matching problem is to identify a set of approximate matches of Q in T (allowing for at most E di erences). The di erences between two sequences of the same species can result from sequencing errors [18, 54] and/or genetic variations [5, 50] . Reads are prone to sequencing errors, which account for about 0.1% of the length of short reads [60, 61, 141] and 10-15% of the length of long reads [19, 83, 165, 170] .
The di erences, known as edits, can be classi ed as substitutions, deletions, or insertions in one or both sequences [100] . Figure 2 shows each possible kind of edit. In ASM, to detect a deleted character or an inserted character, we need to examine all possible pre xes (i.e., substrings that include the rst character of the string) or su xes (i.e., substrings that include the last character of the string) of the two input sequences, and keep track of the pairs of pre xes or su xes that provide the minimum number of edits.
Approximate string matching is needed not only to determine the minimum number of edits between two genomic sequences, but also to provide the location and type of each edit. As two sequences could have a large number of di erent possible arrangements of the edit operations and matches (and hence di erent alignments), the approximate string matching algorithm usually involves a traceback step. The alignment
Read:
insertion substitution deletion C Figure 2 . Three types of errors (i.e., edits). score is the sum of all edit penalties and match scores along the alignment, as de ned by a user-speci ed scoring function. This step nds the optimal alignment as the combination of edit operations to build up the highest alignment score.
Approximate string matching is typically implemented as a dynamic programming based algorithm. Existing implementations, such as Levenshtein distance [100] , Smith-Waterman [154] , and Needleman-Wunsch [126] , have quadratic time and space complexity (i.e., O(m × n) between two sequences with lengths m and n). Therefore, it is desirable to nd lowercomplexity algorithms for ASM.
One candidate to replace dynamic programming based algorithms for ASM is the Bitap algorithm [21, 174] . Bitap tackles the problem of computing the minimum edit distance between a reference text (e.g., reference genome) and a query pattern (e.g., read) with a maximum of k many errors. When k is 0, the algorithm nds the exact matches.
Algorithm 1 shows the Bitap algorithm and Figure 3 shows an example for the execution of the algorithm. The algorithm starts with a pre-processing procedure (Line 4 in Algorithm 1; 0 in Figure 3 ) that converts the query pattern into m-sized pattern bitmasks, PM. We generate one pattern bitmask for each character in the alphabet. Since 0 means match in the Bitap algorithm, we set PM[a][i] = 0 when pattern[i] = a, where a is a character from the alphabet (e.g., A, C, G, T). These pattern bitmasks help us to represent the query pattern in a binary format. After the bitmasks are prepared for each character, every bit of all status bitvectors (R[d], where d is in range [0, k]) is initialized to 1 (Lines 5-6 in Algorithm 1; 0 in Figure 3 ). Each R[d] bitvector at text iteration i holds the partial match information between text[i : (n-1)] (Line 8) and the query with maximum of d errors. Since at the beginning of the execution there are no matches, we initialize all status bitvectors with 1s. The status bitvectors of the previous iteration with edit distance d is kept in oldR[d] (Lines 10-11) to take partial matches into consideration in the next iterations.
The algorithm examines each text character one by one, one per iteration. At each text iteration ( 1 -5 ), the pattern bitmask of the current text character (PM) is retrieved (Line 12). After the status bitvector for exact match is computed (R[0]; Line 13), the status bitvectors for each distance (R[d]; d = 1...k) are computed using the rules in Lines 15-19. For a distance d, three intermediate bitvectors for the error cases (one each for deletion, insertion, substitution; D/I/S in Figure 3 ) are calculated by using oldR [d -1] or R[d -1], since a new error is being added (i.e., the distance is increasing by 1), while the intermediate bitvector for the match case (M) is calculated using oldR [d] . For a deletion (Line 15), we are looking for a string match if the current pattern character is missing, so we copy the partial match information of the previous character (oldR [d -1] ; consuming a text character) without any shifting (not consuming a pattern character) to serve as the deletion bitvector (labeled as D of R1 bitvectors in 1 -5 ) . For a substitution (Line 16), we are looking for a string match if the current pattern character and the current text character do not match, so we take the partial match information of the previous character (oldR [d -1] ; consuming a text character) and shift it left by one (consuming a pattern character) before saving it as the substitution bitvector (labeled as S of R1 bitvectors in 1 -5 text character is missing, so we copy the partial match information of the current character (R[d -1]; not consuming a text character) and shift it left by one (consuming a pattern character) before saving it as the insertion bitvector (labeled as I of R1 bitvectors in 1 -5 ). For a match (Line 18), we are looking for a string match only if the current pattern character matches the current text character, so we take the partial match information of the previous character (oldR [d] ;
consuming a text character but not increasing the edit distance), shift it left by one (consuming a pattern character), and perform an OR operation with the pattern bitmask of the current text character (curPM; comparing the text character and the pattern character) before saving the result as the match bitvector (labeled as R0 bitvectors and M of R1 bitvectors in 1 -5 ).
After computing all four intermediate bitvectors, in order to take all possible partial matches into consideration, we perform an AND operation (Line 19) with these four bitvectors to preserve all 0s that exist in any of them (i.e., all potential locations for a string match with an edit distance of d up to this point). We save the ANDed result as the R[d] status bitvector for the current iteration. This process is repeated for each potential edit distance value from 0 to k. If the most signi cant bit of the R[d] bitvector becomes 0 (Lines 20-22), then there is a match starting at position i of the text with an edit distance d (as shown in 3 -5 ) . The traversal of the text then continues until all possible text positions are examined.
Although the Bitap algorithm is highly suitable for hardware acceleration due to the simple nature of its bitwise operations, we nd that it has ve limitations that hinder its applicability and e cient hardware acceleration for genome analysis. In this section, we discuss each of these limitations.
In order to overcome these limitations and design an e ective and e cient accelerator, we nd that we need to both (1) modify and extend the Bitap algorithm and (2) develop specialized hardware that can exploit the new opportunities that our algorithmic modi cations provide.
No Support for Long Reads. In state-of-the-art implementations of Bitap, the query length is limited by the word size of the machine running the algorithm. This is due to (1) the fact that the bitvector length must be equal to the query length, and (2) the need to perform bitwise operations on the bitvectors. By limiting the bitvector length to a word, each bitwise operation can be done using a single CPU instruction. Unfortunately, the lack of multi-word queries prevents these implementations from working for long reads, whose lengths are on the order of thousands to millions of base pairs (which require thousands of bits to store).
Data . This two-level data dependency forces the consecutive iterations to take place sequentially.
No Support for Traceback. Although the baseline Bitap algorithm can nd possible matching locations of each query read within the reference text, this covers only the rst step of approximate string matching required for genome sequence analysis. Since there could be multiple di erent alignments between the read and the reference, the traceback operation [14, 51, 62, 63, 117, 120, 154, 163, 168, 169] is needed to nd the optimal alignment, which is the alignment with the minimum edit distance (or with the highest score based on a user-de ned scoring function). However, Bitap does not include any such support for optimal alignment identi cation. Limited Compute Parallelism. Even after we solve the algorithmic limitations of Bitap, we nd that we cannot extract signi cant performance bene ts with just algorithmic enhancements alone. For example, while Bitap iterates over each character of the input text sequentially (Line 8), we can enable text-level parallelism to improve its performance (Section 5). However, the achievable level of parallelism is limited by the number of compute units in existing systems. For example, our studies show that Bitap is bottlenecked by computation on CPUs, since the working set ts within the private caches but the limited number of cores prevents the further speedup of the algorithm.
Limited Memory Bandwidth. We would expect that a GPU, which has thousands of compute units, can overcome the limited compute parallelism issues that CPUs experience. However, we nd that a GPU implementation of the Bitap algorithm su ers from the limited amount of memory bandwidth available for each GPU thread. Even when we run a CUDA implementation of the baseline Bitap algorithm [104] , whose bandwidth requirements are signi cantly lower than our modi ed algorithm, the limited memory bandwidth bottlenecks the algorithm's performance. We nd that the bottleneck is exacerbated after the number of threads per block reaches 32, as Bitap becomes shared cache-bound (i.e., on-GPU L2 cache-bound). The small number of registers becomes insu cient to hold the intermediate data required for Bitap execution. Furthermore, when the working set of a thread does not t within the private memory of the thread, destructive interference between threads while accessing the shared memory creates bottlenecks in the algorithm on GPUs. We expect these issues to worsen when we implement traceback, which requires signi cantly higher bandwidth than Bitap.
Our goal in this work is to overcome these limitations and use Bitap in a fast, e cient, and exible ASM framework for both short and long reads. We nd that this goal cannot be achieved by modifying only the algorithm or only the hardware. We design GenASM, the rst ASM acceleration framework for genome sequence analysis. Through careful modi cation and co-design of the enhanced Bitap algorithm and hardware, GenASM aims to successfully replace the expensive dynamic programming based algorithm used for ASM in genomics with the e cient bitwise-operation-based Bitap algorithm, which can accelerate multiple steps of genome sequence analysis.
In GenASM, we co-design our modi ed Bitap algorithm for distance calculation (DC) and our new Bitap-compatible traceback (TB) algorithm with an area-and power-e cient hardware accelerator. GenASM consists of two components, as shown in Figure 4 : (1) GenASM-DC (Section 5), which for each read generates the bitvectors and performs the minimum edit distance calculation (DC); and (2) GenASM-TB (Section 6), which uses the bitvectors to perform traceback (TB) and nd the optimal alignment. GenASM is a exible framework that can be used for di erent use cases (Section 8).
GenASM execution starts when the host CPU issues a task to GenASM with the reference and the query sequences' locations ( 1 in Figure 4 ). GenASM-DC reads the corresponding reference text region and the query pattern from the memory. GenASM-DC then writes these to its dedicated SRAM, which we call DC-SRAM ( 2 ). After that, GenASM-DC divides the reference text (e.g., reference genome) and query pattern (e.g., read) into multiple overlapping windows ( 3 ), and for each sub-text (i.e., the portion of the reference text in one window) and sub-pattern (i.e., the portion of the query pattern in one window), GenASM-DC searches for the sub-pattern within the sub-text and generates the bitvectors ( 4 ). Each processing element (PE) of GenASM-DC writes the generated bitvectors to its own dedicated SRAM, which we call TB-SRAM ( 5 ). Once GenASM-DC completes its search for the current window, GenASM-TB starts reading the stored bitvectors from TB-SRAMs ( 6 ) and generates the window's traceback output ( 7 ). Once GenASM-TB generates this output, GenASM computes the next window and repeats Steps 3 -7 until all windows are completed.
Our hardware accelerators are designed to maximize parallelism and minimize memory footprint. Our modi ed GenASM-DC algorithm is highly parallelizable, and performs only simple and regular bitwise operations, so we implement the GenASM-DC accelerator as a systolic array based accelerator. GenASM-TB accelerator requires simple logic operations to perform the TB-SRAM accesses and the required control ow to complete the traceback operation. Both of our hardware accelerators are highly e cient in terms of area and power. We discuss them in detail in Section 7.
We modify the baseline Bitap algorithm (Section 2.3) to (1) enable e cient alignment of long reads, (2) remove the data dependency between the iterations, and (3) provide parallelism for the large number of iterations.
Long Read Support. The GenASM-DC algorithm overcomes the word-length limit of Bitap (Section 3.1) by storing the bitvectors in multiple words when the query is longer than the word size. Although this modi cation leads to additional computation when performing shifts, it helps GenASM to support both short and long reads. When shifting word i of a multi-word bitvector, the bit shifted out (MSB) of word i -1 needs to be stored separately before performing the shift on word i -1. Then, that saved bit needs to be loaded as the least signi cant bit (LSB) of word i when the shift occurs. This causes the complexity of the algorithm to be m w × n × k, where m is the query length, w is the word size, n is the text length, and k is the edit distance.
Loop Dependency Removal. In order to solve the twolevel data dependency limitation of the baseline Bitap algorithm (Section 3.1), GenASM-DC performs loop unrolling and enables computing non-neighbor (i.e., independent) bitvectors in parallel. Figure 5 shows an example for unrolling with four threads for text characters T0-T3 and status bitvectors R0-R7. For the iteration where R[d] represents T2-R2 (i.e., the target cell shaded in dark red), R[d -1] refers to T2-R1, oldR[d -1] refers to T1-R1, and oldR[d] refers to T1-R2 (i.e., cells T2-R2 is dependent on, shaded in light red). Based on this example, T2-R2 depends on T1-R2, T2-R1, and T1-R1, but it does not depend on T3-R1, T1-R3, or T0-R4. Thus, these independent bitvectors can be computed in parallel without waiting for one another.
Thread1 T2-R0  T1-R1  T0-R2  −  #4  T3-R0  T2-R1  T1-R2  T0-R3  #5  T0-R4  T3-R1  T2-R2  T1-R3  #6  T1-R4  T0-R5  T3-R2  T2-R3  #7  T2-R4  T1-R5  T0-R6  T3-R3  #8  T3-R4  T2-R5  T1-R6  T0-R7  #9  −  T3-R5  T2- 
data written to memory data read from memory
T3-R7 Text-Level Parallelism. In addition to the parallelism enabled by removing the loop dependencies, we enable GenASM-DC algorithm to exploit text-level parallelism. This parallelism is enabled by dividing the text into overlapping sub-texts and searching the query in each of these sub-texts in parallel. The overlap ensures that we do not miss any possible match that may fall around the edges of a sub-text. To guarantee this, the overlap needs to be of length m + k, where m is the query length and k is the edit distance threshold.
After nding the matching location of the text and the edit distance with GenASM-DC, our new traceback [14, 51, 62, 63, 117, 120, 154, 163, 168, 169] algorithm, GenASM-TB, nds the sequence of matches, substitutions, insertions and deletions, along with their positions (i.e., CIGAR string) for the matched region (i.e., the text region that starts from the location reported by GenASM-DC and has a length of m + k), and reports the optimal alignment. Traceback execution (1) starts from the rst character of the matched region between the reference text and query pattern, (2) examines each character and decides which of the four operations should be picked in each iteration, and (3) ends when we reach the last character of the matched region. GenASM-TB uses the intermediate bitvectors generated and saved in each iteration of the GenASM-DC algorithm (i.e., match, substitution, deletion and insertion bitvectors generated in Lines 15-18 in Algorithm 1). After a value 0 is found at the MSB of one of the R[d] bitvectors (i.e., a string match is found with d errors), GenASM-TB walks through the bitvectors back to the LSB, following a chain of 0s (which indicate matches at each location) and reverting the bitwise operations. At each position, based on which of the four bitvectors holds a value 0 in each iteration (starting with an MSB with a 0 and ending with an LSB with a 0), the sequence of matches, substitutions, insertions and deletions (i.e., traceback output) is found for each position of the corresponding alignment found by GenASM-DC. Unlike GenASM-DC, GenASM-TB has an irregular control ow within the stored intermediate bitvectors, which depends on the text and the pattern.
Algorithm 2 shows the GenASM-TB algorithm and Figure 6 shows an example for the execution of the algorithm for each of the alignments found in 3 -5 of Figure 3 . In Figure 6 , <x, y, z> stands for patternI, textI and curError, respectively (Lines 6-8 in Algorithm 2). patternI represents the position of a 0 currently being processed within a given bitvector (i.e., pattern index), textI represents the outer loop iteration index (i.e., text index; i in Algorithm 1), and curError represents the inner loop iteration index (i.e., number of remaining errors; d in Algorithm 1).
When we nd a 0 at match[textI][curError][patternI] (i.e., a match (M) is found for the current position; Line 17), one character each from both text and query is consumed, but the number of remaining errors stays the same. Thus, the pointer moves to the next text character (as the text character is consumed), and the 0 currently being processed (highlighted with orange color in Figure 6 ) is right-shifted by one (as the query character is also consumed). In other words, textI is incremented (Line 28), patternI is decremented (Line 30), but curError remains the same. Thus, <x, y, z> becomes <x -1, y + 1, z> after we nd a match. For example, in Figure 6a , for Text[0], we have <3, 0, 1> for the indices, and after the match is found, at the next position (Text[1]), we have <2, 1, 1>.
When we nd a 0 at subs[textI][curError][patternI] (i.e., a substitution (S) is found for the current position; Line 19), one character each from both text and query is consumed, and the number of remaining errors is decremented (Line 26). Thus, <x, y, z> becomes <x -1, y + 1, z -1> after we nd a substitution (e.g., Text [1] in Figure 6b) .
When we nd a 0 at ins[textI][curError][patternI] (i.e., an insertion (I) is found for the current position; Lines 13 and 21), the inserted character does not appear in the text, and only a character from the pattern is consumed. The 0 currently being processed is right-shifted by one, but the text pointer remains the same, and the number of remaining errors is decremented. Thus, <x, y, z> becomes <x -1, y, z -1> after we nd an insertion (e.g., Text [-] in Figure 6c) .
When we nd a 0 at del[textI][curError][patternI] (i.e., a deletion (D) is found for the current position; Lines 15 and 23), the deleted character does not appear in the pattern, and only a character from the text is consumed. The 0 currently being processed is not right-shifted, but the pointer moves to the next text character, and the number of remaining errors is also decremented. Thus, <x, y, z> becomes <x, y + 1, z -1> after we nd an insertion (e.g., Text [1] in Figure 6a ). Divide-and-Conquer Approach. Since GenASM-DC stores all of the intermediate bitvectors, in the worst case, the length of the text region that the query pattern maps to can be m + k, assuming all of the errors are deletions from the pattern. Since we need to store all of the bitvectors for m + k characters, and compute 4 × k many bitvectors within each text iteration (each m bits long), for long reads with high error rates, the memory requirement becomes~80GB, when m is 10,000 and k is 1,500.
In order to decrease the memory footprint of the algorithm, we follow two key ideas. First, we apply a divide-and-conquer approach (similar to the tiling approach of Darwin's alignment accelerator, GACT [162] ). Instead of storing all of the bitvectors for m + k text characters, we divide the text and pattern into overlapping windows (i.e., sub-text and sub-pattern; Lines 3-4 in Algorithm 2) and perform the traceback computation for each window. After all of the windows' partial traceback outputs are generated, we merge them to nd the complete traceback output. This approach helps us to decrease the memory footprint from ((m + k) × 4 × k × m) bits to (W × 4 × W × W ) bits, where W is the window size. This divide-and-conquer approach also helps us to reduce the complexity of the bitvector generation step (Section 5) from m w × n × k to W w × W × W . Second, instead of storing all 4 bitvectors (i.e., match, substitution, insertion, deletion) separately, we only need to store bitvectors for match, insertion, and deletion, as the substitution bitvector can be obtained easily by left-shifting the deletion bitvector by 1 (Line 16 in Algorithm 1). This modi cation helps us to decrease the required write bandwidth and the memory footprint to (W × 3 × W × W ) bits.
GenASM-TB restricts the number of consumed characters from the text or the pattern to W-O (Line 11 in Algorithm 2) to ensure that consecutive windows share O characters (i.e., overlap size), and thus, the traceback output can be generated accurately. The sub-text and the sub-pattern corresponding to each window are found using the number of consumed text characters (textConsumed) and the number of consumed pattern characters (patternConsumed) in the previous window (Lines 31-32 in Algorithm 2).
Partial Support for Complex Scoring Schemes. We extend the GenASM-TB algorithm to provide partial support (Section 10.2) for non-unit costs for di erent edits and the a ne gap penalty model [14, 62, 117, 168] . By changing the order in which di erent traceback cases are checked in Lines 13-24 in Algorithm 2, we can support di erent types of scoring schemes. For example, in order to mimic the behavior of the a ne gap penalty model, we check whether the traceback output that has been chosen for the previous position (i.e., prev) is an insertion or a deletion. If the previous edit is a gap (insertion or deletion), and there is a 0 at the current position of the insertion or deletion bitvector (Lines 13 and 15 in Algorithm 2), then we prioritize extending this previously opened gap, and choose insertion-extend or deletion-extend as the current position's traceback output, depending on the type of the previous gap. As another example, in order to mimic the behavior of non-unit costs for di erent edits, we can simply sort three error cases 
GenASM-DC Hardware. We implement GenASM-DC as a linear cyclic systolic array [93, 94] based accelerator. The accelerator is optimized to reduce both the memory bandwidth and the memory footprint. Feedback logic enabling cyclic systolic behavior allows us to x the required number of memory ports [93] and to reduce memory footprint.
A GenASM-DC accelerator consists of a processing block (PB; Figure 7a ) along with a control and memory management logic. A PB consists of multiple processing elements (PEs). Each PE contains a single processing core (PC; Figure 7b ) and ip-op-based storage logic. The PC is the primary compute unit, and implements Lines 15-19 of Algorithm 1 to perform the approximate string matching for a w-bit query pattern. The number of PEs in a PB is based on compute, area, memory bandwidth and power requirements. This block also implements the logic to load data from outside of the array (i.e., DC-SRAM; Figure 7a ) or internally for cyclic operations.
GenASM-DC uses two types of SRAM bu ers (Figure 7a ): (1) DC-SRAM, which stores the reference text, the pattern bitmasks for the query read, and the intermediate data generated from PEs (i.e., oldR values and MSBs required for shifts; Section 5); and (2) TB-SRAM, which stores the intermediate bitvectors from GenASM-DC for later use by GenASM-TB. For a 64-PE con guration with 64 bits of processing per PE, and for the case where we have a long (10Kbp) read 1 with a high error rate (15%) and a corresponding text region of 11.5Kbp, GenASM-DC requires a total of 8KB DC-SRAM storage. For each PE, we have a dedicated TB-SRAM, which stores the match, insertion and deletion bitvectors generated by the corresponding PE. For the same con guration of GenASM-DC, each PE requires a total of 1.5KB TB-SRAM storage, with a single R/W port. In each cycle, 192 bits of data (24B) is written to each TB-SRAM by each PE.
When each thread (i.e., each column) in Figure 5 is mapped to a PE, GenASM-DC coordinates the data dependencies across DC iterations, with the help of two ip-ops in each PE. For example, T2-R2 in Figure 5 is generated by PE x in Cycle y , and is mapped to R[d]. In order to generate T2-R2, T2-R1 (which maps to R[d -1]) needs to be generated by PE x-1 in Cycle y-1 ( 1 in Figure 7 ), T1-R1 (which maps to oldR[d -1]) needs to be generated by PE x-1 in Cycle y-2 ( 2 ), and T1-R2 (which maps to oldR[d]) needs to be generated by PE x in Cycle y-1 ( 3 ), where x is the PE index and y is the cycle index. With this dependency-aware mapping, regardless of the number of instantiated PEs, we can successfully limit DC-SRAM tra c for a single PB to only one read and one write per cycle.
GenASM-TB Hardware. After GenASM-DC nishes writing all of the intermediate bitvectors to TB-SRAMs, GenASM-TB reads them by following an irregular control ow, which depends on the text and the pattern to nd the optimal alignment (by implementing Algorithm 2).
In our GenASM con guration, where we have 64 PEs and 64 bits per PE in a GenASM-DC accelerator, and the window size (W ) is 64 (Section 6), we have one 1.5KB TB-SRAM (which ts our 24B/cycle × 64 cycles/window output storage requirement) for each of the 64 PEs. As Figure 8 shows, a single GenASM-TB accelerator is connected to all of these 64 TB-SRAMs (96KB, in total). In each GenASM-TB cycle, we read from only one TB-SRAM. curError provides the 1 Although we use 10Kbp-long reads in our analysis (Section 9), GenASM does not have any limitation on the length of reads as a result of our divideand-conquer approach (Section 6). index of the TB-SRAM that we read from; textI provides the starting index within this TB-SRAM, which we read the next set of bitvectors from; and patternI provides the position of the 0 being processed (Algorithm 2). We implement the GenASM-TB hardware using very simple logic (Figure 8) , which 1 reads the bitvectors from one of the TB-SRAMs using the computed address, 2 performs the required bitwise comparisons to nd the CIGAR character for the current position, and 3 computes the next TB-SRAM address to read the new set of bitvectors. After GenASM-TB nds the complete CIGAR string, it writes the output to main memory and completes its execution.
PE2 PE65 Overall System. We design our system to take advantage of modern 3D-stacked memory systems [58, 92] , such as the Hybrid Memory Cube (HMC) [76] or High-Bandwidth Memory (HBM) [86, 99] . Such memories are made up of multiple layers of DRAM arrays that are stacked vertically in a single package. These layers are connected via high-bandwidth links called through-silicon vias (TSVs) that provide lower-latency and more energy-e cient data access to the layers than the external DRAM I/O pins [39, 99] . Memories such as HMC and HBM include a dedicated logic layer that connects to the TSVs and allows processing elements to be implemented in memory to exploit the e cient data access. Due to thermal and area constraints, only simple processing elements that execute low-complexity operations (e.g., bitwise logic, simple arithmetic, simple cores) can be included in the logic layer [3, 4, 23, 24, 43, 56, 72, 73, 91, 119, 137] .
We decide to implement GenASM in the logic layer of 3Dstacked memory, for two reasons. First, we can exploit the natural subdivision within 3D-stacked memory (e.g., vaults in HMC [76] , pseudo-channels in HBM [86] ) to e ciently enable parallelism across multiple GenASM accelerators. This subdivision allows accelerators to work in parallel without interfering with each other. Second, we can reduce the power consumed for DRAM accesses by reducing o -chip data movement across the memory channel [119] . Both of our hardware accelerators are highly e cient in terms of area and power (Section 10.1), and can t within the logic layer's constraints.
To illustrate how GenASM takes advantage of 3D-stacked memory, we discuss an example implementation of GenASM inside the logic layer of a 16GB HMC with 32 vaults [76] . Within each vault, the logic layer contains a GenASM-DC accelerator, its associated DC-SRAM (8KB), a GenASM-TB accelerator, and TB-SRAMs (64×1.5KB). Since we have small SRAM bu ers for both DC and TB to exploit locality, GenASM accesses the memory and utilizes the memory bandwidth only to read the reference and the query sequences. One GenASM accelerator at each vault requires 105-142 MB/s bandwidth, thus the total bandwidth requirement of all 32 GenASM accelerators is 3.3-4.4 GB/s (which is much less than peak bandwidth provided by modern 3D-stacked memories).
We demonstrate the e ciency and exibility of the GenASM acceleration framework by describing three use cases of approximate string matching in genome sequence analysis: (1) read alignment step of short and long read mapping, (2) pre-alignment ltering for short reads, and (3) edit distance calculation between any two sequences. We believe the GenASM framework can be useful for many other use cases, and we discuss some of them brie y in Section 11.
Read Alignment of Short and Long Reads. As we explain in Section 2.1, read alignment is the last step of short and long read mapping. In read alignment, all of the remaining candidate mapping regions of the reference genome and the query reads are aligned, in order to identify the mapping that yields either the lowest total number of errors (if using edit distance based scoring) or the highest score (if using a user-de ned scoring function). Thus, read alignment can be a use case for approximate string matching, since errors (i.e., substitutions, insertions, deletions) should be taken into account when aligning the sequences. As part of read alignment, we also need to generate the traceback output for the best alignment between the reference region and the read.
For read alignment, the whole GenASM pipeline, as explained in Section 4, should be executed, including the traceback step. In general, read alignment requires more complex scoring schemes, where di erent types of edits have non-unit costs. Thus, GenASM-TB should be con gured based on the given cost of each type of edit (Section 6). As GenASM framework can work with arbitrary length sequences, we can use it to accelerate both short read and long read alignment.
Pre-Alignment Filtering for Short Reads. In the prealignment ltering step of short read mapping, the candidate mapping locations, reported by the seeding step, are further ltered by using di erent mechanisms. Although the regions of the reference at these candidate mapping locations share common seeds with query reads, they are not necessarily similar sequences. To avoid examining dissimilar sequences at the downstream computationally-expensive read alignment step, a pre-alignment lter estimates the edit distance between every read and the regions of the reference at each read's candidate mapping locations, and uses this estimation to quickly decide whether or not read alignment is needed. If the sequences are dissimilar enough, signi cant amount of time is saved by avoiding the expensive alignment step [9, 10, 13, 176, 177] .
In pre-alignment ltering, since we only need to estimate (approximately) the edit distance and check whether it is above a user-de ned threshold, GenASM-DC can be used as a pre-alignment lter. As GenASM-DC is very e cient when we have shorter sequences and a low error threshold (due to the O(m × n × k) complexity of the underlying Bitap algorithm, where m is the query length, n is the reference length, and k is the number of allowed errors), GenASM framework can e ciently accelerate the pre-alignment ltering step of especially short read mapping. 2 Edit Distance Calculation. Edit distance, also called Levenshtein distance [100] , is the minimum number of edits (i.e., substitutions, insertions and deletions) required to convert one sequence to another. Edit distance calculation is one of the fundamental operations in genomics to measure the similarity or distance between two sequences [155] . As we explain in Section 2.3, the Bitap algorithm, which is the underlying algorithm of GenASM-DC, is originally designed for edit distance calculation. Thus, GenASM framework can accelerate edit distance calculation between any two arbitrary-length genomic sequences.
Although GenASM-DC can nd the edit distance by itself and traceback is optional for this use case, DC-TB interaction is required in our accelerator to exploit the e cient divideand-conquer approach GenASM follows. Thus, GenASM-DC and GenASM-TB work together to nd the minimum edit distance in a fast and memory-e cient way, but the traceback output is not generated or reported by default (though it can optionally be enabled).
Area and Power Analysis. We synthesize and place & route the GenASM-DC and GenASM-TB accelerator datapaths using the Synopsys Design Compiler [156] with a typical 28nm low-power process, with memories generated using an industry-grade SRAM compiler, to analyze the accelerators' area and power. Our synthesis targets post-routing timing closure at 1GHz clock frequency. We then use an in-house cycle-accurate simulator parameterized with the synthesis and memory estimations to drive the performance and power analysis.
We evaluate a 16GB HMC-like 3D-stacked DRAM architecture, with 32 vaults [76] and 256GB/s of internal bandwidth [23, 76] , and a clock frequency of 1.25GHz [76] . The amount of available area in the logic layer for GenASM is around 3.5-4.4 mm 2 per vault [23, 43] . The power budget of our PIM logic per vault is 312mW [43] .
Performance Model. We build a spreadsheet-based analytical model for GenASM-DC and GenASM-TB, which considers reference genome (i.e., text) length, query read (i.e., pattern) length, maximum edit distance, window size, hardware design parameters (number of PEs, bit width of each PE) and number of vaults as input parameters and projects compute cycles, DRAM read/write bandwidth, SRAM read/write bandwidth, and memory footprint. We verify the analyticallyestimated cycle counts for various PE con gurations with the cycle counts collected from our RTL simulations.
Read Alignment Comparisons. For the read alignment use case, we compare GenASM with the read alignment steps of two commonly-used state-of-the-art read mappers: Min-imap2 [102] and BWA-MEM [101] , running on an Intel ® Xeon ® Gold 6126 CPU [80] operating at 2.60GHz, with 64GB DDR4 memory. Software baselines are run with a single thread and with 12 threads. We measure the execution time and power consumption of the alignment steps in Minimap2 and BWA-MEM. We measure the individual power consumed by each tool using Intel's PCM power utility [81] .
We also compare GenASM with a state-of-the-art GPUaccelerated short read alignment tool, GASAL2 [2] . We run GASAL2 on an Nvidia Titan V GPU [129] with 12GB HBM2 memory [86] . To fully utilize the GPU, we con gure the number of alignments per batch based on the GPU's number of multiprocessors and the maximum number of threads per multiprocessor, as described in the GASAL2 paper [2] . To better analyze the high parallelism that the GPU provides, we replicate our datasets to obtain datasets with 100K, 1M and 10M reference-read pairs for short reads. We run the datasets with GASAL2, and collect kernel time and average power consumption using nvprof [130] .
We also compare GenASM with two state-of-the-art hardware-based alignment accelerators, GACT of Darwin [162] and SillaX of GenAx [55] . We synthesize and execute the open-source RTL for GACT [161] . We estimate the performance of SillaX using data reported by the original work [55] .
We analyze the alignment accuracy of GenASM by comparing the alignment outputs (i.e., alignment score, edit distance, and CIGAR string) of GenASM with the alignment outputs of BWA-MEM and Minimap2, for short reads and long reads, respectively. We obtain the BWA-MEM and Minimap2 alignments by running the tools with their default settings.
Pre-Alignment Filtering Comparisons. We compare GenASM with Shouji [9] , which is the state-of-the-art FPGAbased pre-alignment lter for short reads. For execution time and ltering accuracy analyses, we use data reported by the original work [9] . For power analysis, we report the total power consumption of Shouji using the power analysis tool in Xilinx Vivado [175] , after synthesizing and implementing the open-source FPGA design of Shouji [149] .
Edit Distance Calculation Comparisons. We compare GenASM with the state-of-the-art software-based read alignment library, Edlib [155] , running on an Intel ® Xeon ® Gold 6126 CPU [80] operating at 2.60GHz, with 64GB DDR4 memory. Edlib uses the Myers' bitvector algorithm [121] to nd the edit distance between two sequences. We use the default global Needleman-Wunsch (NW) [126] mode of Edlib to perform our comparisons. We measure the power consumed by Edlib using Intel's PCM power utility [81] .
We also compare GenASM with ASAP [22] , which is the state-of-the-art FPGA-based accelerator for computing the edit distance between two short reads. We estimate the performance of ASAP using data reported by the original work [22] .
Datasets. For the read alignment use case, we evaluate GenASM using the latest major release of the human genome assembly, GRCh38 [124]. We use the 1-22, X, and Y chromosomes by ltering the unmapped contigs, unlocalized contigs, and mitochondrial genome. Genome characters are encoded into 2-bit patterns (A = 00, C = 01, G = 10, T = 11). With this encoding, the reference genome uses 715 MB of memory.
We generate four sets of long reads (i.e., PacBio and ONT datasets) using PBSIM [131] and three sets of short reads (i.e., Illumina datasets) using Mason [71] . For the PacBio datasets, we use the default error pro le for the continuous long reads (CLR) in PBSIM. For the ONT datasets, we modify the settings to match the error pro le of ONT reads sequenced using R9.0 chemistry [84] . Both datasets have 240,000 reads of length 10Kbp, each simulated with 10% and 15% error rates. The Illumina datasets have 200,000 reads of length 100bp, 150bp, and 250bp, each simulated with a 5% error rate.
For the pre-alignment ltering use case, we use two datasets that Shouji [9] provides as test cases: reference-read pairs (1) of length 100bp with an edit distance threshold of 5, and (2) of length 250bp with an edit distance threshold of 15.
For the edit distance calculation use case, we use the publicly-available dataset that Edlib [155] provides. The dataset includes two real DNA sequences, which are 100Kbp and 1Mbp in length, and arti cially-mutated versions of the original DNA sequences with measures of similarity ranging between 60%-99%. Evaluating this set of sequences with varying values of similarity and length enables us to demonstrate how these parameters a ect performance. Table 1 shows the area and power breakdown of each component in GenASM, and the total area overhead and power consumption of (1) a single GenASM accelerator (in 1 vault) and (2) 32 GenASM accelerators (in 32 vaults). Both GenASM-DC and GenASM-TB operate at 1GHz.
The area overhead of one GenASM accelerator is 0.334 mm 2 , and the power consumption of one GenASM accelerator, including the SRAM power, is 101 mW. When we compare GenASM with a single core of a modern Intel ® Xeon ® Gold 6126 CPU [80] (which we conservatively estimate to use 10.4 W [80] and 32.2 mm 2 [36] per core), we nd that GenASM is signi cantly more e cient in terms of both area and power consumption. As we have one GenASM accelerator per vault, the total area overhead of GenASM in all 32 vaults is 10.69 mm 2 . Similarly, the total power consumption of 32 GenASM accelerators is 3.23 W. 
Software Baselines (CPU). Figure 9 shows the read alignment throughput (reads/sec) of GenASM and the alignment steps of BWA-MEM and Minimap2, when aligning long noisy PacBio and ONT reads against the human reference genome. When comparing with BWA-MEM, we run GenASM with the candidate locations reported by BWA-MEM's ltering step. Similarly, when comparing with Minimap2, we run GenASM with the candidate locations reported by Minimap2's ltering step. GenASM's throughput is determined by the throughput of the execution of GenASM-DC and GenASM-TB with window size (W ) of 64 and overlap size (O) of 24.
As Figure 9 shows, GenASM provides (1) 7173× and 648× throughput improvement over the alignment step of BWA-MEM for its single-thread and 12-thread execution, respectively, and (2) 1126× and 116× throughput improvement over the alignment step of Minimap2 for its single-thread and 12-thread execution, respectively. Based on our power analysis with long reads, we nd that power consumption of BWA-MEM's alignment step is 58.6 W and 109.5 W, and power consumption of Minimap2's read alignment step is 59.8 W and 118.9 W for their single-thread and 12-thread executions, respectively. GenASM consumes only 3.23W, and thus reduces the power consumption of the alignment steps of BWA-MEM and Minimap2 by 18× and 19× for single-thread execution, and by 34× and 37× for 12-thread execution, respectively. Figure 10 compares the read alignment throughput (reads/sec) of GenASM with that of the alignment steps of BWA-MEM and Minimap2, when aligning short Illumina reads against the human reference genome. GenASM provides (1) 1390× and 111× throughput improvement over the alignment step of BWA-MEM for its single-thread and 12-thread execution, respectively, and (2) 1839× and 158× throughput improvement over the alignment step of Min-imap2 for its single-thread and 12-thread execution. Based on our power analysis with short reads, we nd that GenASM reduces the power consumption over the alignment steps of BWA-MEM and Minimap2 by 16× and 18× for single-thread execution, and by 33× and 31× for 12-thread execution, respectively. Figure 11 shows the total execution time of the entire BWA-MEM and Minimap2 pipelines, along with the total execution time when the alignment steps of each pipeline are replaced by GenASM, for the three representative input datasets. As Figure 11 shows, GenASM provides (1) 2.4× and 1.9× speedup for Illumina reads (250bp); (2) 6.5× and 3.4× speedup for PacBio reads (15%); and (3) 4.9× and 2.1× speedup for ONT reads (15%), over the entire pipeline executions of BWA-MEM and Minimap2, respectively. Software Baselines (GPU). We compare GenASM with the state-of-the-art GPU aligner, GASAL2 [2] , using three datasets of varying size (100K, 1M, and 10M reference-read pairs). Based on our analysis, we make three ndings. First, for 100bp Illumina reads, GenASM provides 9.9×, 9.2×, and 8.5× speedup over GASAL2, while reducing the power consumption by 15 20 .2×, and 20.6× for 100K, 1M, and 10M datasets, respectively. We conclude that GenASM provides signi cant performance bene ts and energy e ciency over GPU aligners for short reads.
Hardware Baselines. We compare GenASM with two state-of-the-art hardware accelerators for read alignment: GACT (from Darwin [162] ) and SillaX (from GenAx [55] ).
Darwin is a hardware accelerator designed for long read alignment [162] . Darwin contains components that accelerate both the ltering (D-SOFT) and alignment (GACT) steps of read mapping. The open-source RTL code available for the GACT accelerator of Darwin allows us to estimate the throughput, area and power consumption of GACT and compare it with GenASM for read alignment. In Darwin, GACT logic and the associated 128KB SRAM are responsible for lling the dynamic programming matrix, generating the traceback pointers and nding the maximum score. Thus, we believe that it is fair to compare the power consumption and the area of the GACT logic and GenASM logic, along with their associated SRAMs.
In order to have an iso-bandwidth comparison with Darwin's GACT, we compare only a single array of GACT and a single set of GenASM-DC and GenASM-TB, because (1) GenASM utilizes the high memory bandwidth that PIM provides only to parallelize many sets of GenASM-DC and GenASM-TB, and a single set of GenASM-DC and GenASM-TB does not require high bandwidth, and (2) all internal data of both GenASM and Darwin is provided by local SRAMs. We synthesize both designs (i.e., GenASM and GACT) at an iso-PVT (process, voltage, temperature) corner, with the same number of PEs, and with their optimum parameters.
As Figure 12 shows, for a single GACT array with 64 PEs at 1GHz, the throughput of GACT decreases from 55,556 to 6,289 alignments per second when the sequence length increases from 1Kbp to 10Kbp, while consuming 277.7 mW of power. In comparison, for a single GenASM accelerator at 1GHz (with a 64-PE con guration), the throughput decreases from 236,686 to 23,669 alignments per second when the sequence length increases from 1Kbp to 10Kbp, while consuming 101 mW of power. This shows that, on average, GenASM provides 3.9× better throughput than GACT, while consuming 2.7× less power. Thus, GenASM provides 10.5× better throughput per unit power for long reads when compared to GACT. As Figure 13 shows, we also compare the throughput of GenASM and GACT for short read alignment (i.e., 100-300bp reads). We nd that GenASM performs 7.4× better than GACT when aligning short reads, on average. Thus, GenASM provides 20.0× better throughput per unit power for short reads when compared to GACT. We compare the required area for the GACT logic with 128KB of SRAM and the required area for the GenASM logic (GenASM-DC and GenASM-TB) with 8KB of DC-SRAM and 96KB of TB-SRAMs, at 28nm. We nd that GenASM requires 1.7× less area than GACT. Thus, GenASM provides 6.6× and 12.6× better throughput per unit area for long reads and for short reads, respectively, when compared to GACT.
The main di erence between GenASM and GACT is the underlying algorithms. GenASM uses our modi ed Bitap algorithm, which requires only simple and fast bitwise operations. On the other hand, GACT uses the complex and computationally more expensive dynamic programming based algorithm for alignment. This is the main reason why GenASM is more e cient than GACT of Darwin.
GenAx is a hardware accelerator designed for short read alignment [55] . Similar to Darwin, GenAx accelerates both the ltering and alignment steps of read mapping. Unlike GenAx, whose design is optimized only for short reads, GenASM is more robust and works with both short and long reads. While we are unable to reimplement GenAx, the throughput analysis of SillaX (the alignment accelerator of GenAx) provided by the original work [55] allows us to provide a performance comparison between GenASM and SillaX for short read alignment.
We compare SillaX with GenASM at their optimal operating frequencies (2GHz for SillaX, 1GHz for GenASM), and nd that GenASM provides 1.9× higher throughput for short reads (101bp) than SillaX (whose approximate throughput is 50M alignments per second). Using the area and power numbers reported for the computation logic of SillaX, we nd that GenASM requires 63% less logic area (2.08 mm 2 vs. 5.64 mm 2 ) and 82% less logic power (1.18 W vs. 6.6 W).
In order to compare the total area of SillaX and GenASM, we perform a CACTI-based analysis [172] for the SillaX SRAM (2.02 MB). We nd that the SillaX SRAM consumes an area of 3.47 mm 2 , resulting in a total area of 9.11 mm 2 for Sil-laX. Although GenASM (10.69 mm 2 ) requires 17% more total area than SillaX, we nd that GenASM provides 1.6× better throughput per unit area for short reads than SillaX.
Accuracy Analysis. We compare the traceback outputs of GenASM and (1) BWA-MEM for short reads, (2) Minimap2 for long reads, to assess the accuracy and correctness of GenASM-TB. We nd that the optimum (W , O) setting (i.e., window size and overlap size) for the GenASM-TB algorithm in terms of performance and accuracy is W = 64 and O = 24. With this setting, GenASM completes the alignment of all reads in each dataset, and increasing the window size does not change the alignment output.
For short reads, we use the default scoring setting of BWA-MEM (i.e., match=+1, substitution=-4, gap opening=-6, and gap extension=-1). For 96.6% of the short reads, GenASM nds an alignment whose score is equal to the score of the alignment reported by BWA-MEM. This fraction increases to 99.7% when we consider scores that are within ±4.5% of the scores reported by BWA-MEM.
For long reads, we use the default scoring setting of Min-imap2 (i.e., match=+2, substitution=-4, gap opening=-4, and gap extension=-2). For 99.6% of the long reads with a 10% error rate, GenASM nds an alignment whose score is within ±0.4% of the score of the alignment reported by Minimap2. For 99.7% of the long reads with a 15% error rate, GenASM nds an alignment whose score is within ±0.7% of the score of the alignment reported by Minimap2.
There are two reasons for the di erence between the alignment scores reported by GenASM and the scores reported by the baseline tools. First, GenASM performs traceback for the alignment with the minimum edit distance. However, the baseline can report an alignment that has a higher number of edits but a lower score than the alignment reported by GenASM, when more complex scoring schemes are used. Second, during the TB stage, GenASM follows a xed order at each iteration when picking between substitutions, insertions, or deletions (based on the penalty of each error type). While we pick the error type with the lowest possible cost at the current iteration, another error type with a higher initial cost may lead to a better (i.e., lower-cost) alignment in later iterations, which cannot be known beforehand. 3 Although GenASM is optimized for unit-cost based scoring (i.e., edit distance) and currently provides only partial support for more complex scoring schemes, we show that GenASM framework can still serve as a fast, memory-and powere cient, and quite accurate alternative for read alignment.
We compare GenASM with the state-of-the-art FPGAbased pre-alignment lter for short reads, Shouji [9] , using two datasets provided in [9] . When we compare Shouji (with maximum ltering units) and GenASM for the dataset with 100bp sequences, we nd that GenASM provides 3.7× speedup over Shouji, while reducing power consumption by 1.7×. When we perform the same analysis with 250bp sequences, we nd that GenASM does not provide speedup over Shouji, but reduces power consumption by 1.6×.
In pre-alignment ltering for short reads, only GenASM-DC is executed (Section 8). The complexity of GenASM-DC is O(n × m × k) whereas the complexity of Shouji is O(m × k), where n is the text length, m is the read length, and k is the edit distance threshold. Going from the 100bp dataset to the 250bp dataset, all these three parameters increase linearly. Thus, the speedup of GenASM over Shouji for pre-alignment ltering decreases for datasets with longer reads.
To analyze ltering accuracy, we use Edlib [155] to generate the ground truth edit distance value for each sequence pair in the datasets (similar to Shouji). We evaluate the accuracy of GenASM as a pre-alignment lter by computing its false accept rate and false reject rate (as de ned in [9] ).
The false accept rate [9] is the ratio of the number of dissimilar sequences that are falsely accepted by the lter (as similar) and the total number of dissimilar sequences that are rejected by the ground truth. The goal is to minimize the false accept rate to maximize the number of dissimilar sequences that are eliminated by the lter. For the 100bp dataset with an edit distance threshold of 5, Shouji has a 4% false accept rate, whereas GenASM has a false accept rate of only 0.02%. For the 250bp dataset with an edit distance threshold of 15, Shouji has a 17% false accept rate, whereas GenASM has a false accept rate of only 0.002%. Thus, GenASM provides a very low rate of falsely-accepted dissimilar sequences, and signi cantly improves the accuracy of pre-alignment ltering compared to Shouji.
While Shouji approximates the edit distance, GenASM calculates the actual distance. Although calculation requires more computation than approximation, a computed distance results in a near-zero (0.002%) false accept rate. 4 Thus, GenASM lters more false-positive locations out, leaving fewer candidate locations for the expensive alignment step to process. This greatly reduces the combined execution time of ltering and alignment. Thus, even though GenASM does not provide any speedup over Shouji when ltering the 250bp sequences, its lower false accept rate makes it a better option for this step of the pipeline with greater overall bene ts.
The false reject rate [9] is the ratio of the number of similar sequences that are rejected by the lter (as dissimilar) and the total number of similar sequences that are accepted by the ground truth. The false reject rate should always be equal to 0%. We observe that GenASM always provides a 0% false reject rate, and thus does not lter out similar sequence pairs, as does Shouji.
We compare GenASM with the state-of-the-art edit distance calculation library, Edlib [155] . Figure 14 compares the execution time of Edlib (with and without nding the traceback output) and GenASM when nding the edit distance between two sequences of length 100Kbp, and also two sequences of length 1Mbp, which have similarity ranging from 60% to 99% (Section 9). Since Edlib is a single-thread edit distance calculation tool, for a fair comparison, we compare the throughput of only one GenASM accelerator (i.e., in one vault) with a single-thread execution of the Edlib tool.
As Figure 14 shows, when performing edit distance calculation between two 100Kbp sequences, GenASM provides 22-716× and 146-1458× speedup over Edlib execution without and with traceback, respectively. GenASM has the same execution time for both of the cases. When the sequence length increases from 100Kbp to 1Mbp, the execution time of GenASM increases linearly (since W is constant, but m + k increases linearly). However, due to its quadratic complexity, Edlib cannot scale linearly. Thus, for the edit distance calculation of 1Mbp sequences, GenASM provides 262-5413× and 627-12501× speedup over Edlib execution without and with traceback, respectively.
Although both the GenASM algorithm and Edlib's underlying Myers' algorithm [121] use bitwise operations only for edit distance calculation and exploit bit-level parallelism, the main advantages of the GenASM algorithm come from (1) the divide-and-conquer approach we follow for e cient support for longer sequences, and (2) our e cient co-design of the GenASM algorithm with the GenASM hardware accelerator. Based on our power analysis, we nd that power consumption of Edlib is 55.3 W and 58.8 W when nding the edit distance between two 100Kbp sequences and two 1Mbp sequences, respectively. Thus, GenASM reduces power consumption by 548× and 582× over Edlib, respectively.
We also compare GenASM with ASAP [22] , the state-ofthe-art FPGA-based accelerator for edit distance calculation. While we are unable to reimplement ASAP, the execution time and power consumption analysis of ASAP provided in [22] allows us to provide a comparison between GenASM and ASAP. ASAP is optimized for shorter sequences and reports execution time only for sequences of length 64bp-320bp [22] . Based on [22] , the execution time of one ASAP accelerator increases from 6.8 µs to 18.8 µs when the sequence length increases from 64bp to 320bp, while consuming 6.8 W of power. In comparison, we report that the execution time of one GenASM accelerator increases from 0.017 µs to 2.025 µs when the sequence length increases from 64bp to 320bp, while consuming 0.101 W of power. This shows that GenASM provides 9.3-400× speedup over ASAP, while consuming 67× less power.
GenASM's performance improvements come from our algorithm/hardware co-design, i.e., both from our modi ed algorithm and our co-designed architecture for this algorithm. The sources of the large improvements in GenASM are (1) the very simple computations it performs; (2) the divideand-conquer approach we follow, which makes our design e cient for both short and long reads despite their di erent error pro les; and (3) the very high degree of parallelism obtained with the help of specialized compute units, dedicated SRAMs for both GenASM-DC and GenASM-TB, and the vault-level parallelism provided by processing in the logic layer of 3D-stacked memory.
Algorithm-Level. Our divide-and-conquer approach allows us to decrease the execution time of GenASM-DC from ( m×(m+k)×k P×w ) cycles to (( W ×W ×min(W ,k) P×w )× m+k W -O ) cycles, where m is the pattern size, k is the edit distance threshold, P is the number of PEs that GenASM-DC has (i.e., 64) , w is the number of bits processed by each PE (i.e., 64) , W is the window size (i.e., 64) , and O is the overlap size between windows (i.e., 24) . Although the total GenASM-TB execution time does not change ((m + k) cycles vs. ((W -O) × m+k W -O ) cycles), our divide-and-conquer approach helps us decrease the GenASM-DC execution time by 3662× for long reads, and by 1.6 -3.9× for short reads.
Hardware-Level. GenASM-DC's systolic-array-based design removes the data dependency limitation of the underlying Bitap algorithm, and provides 64× parallelism by performing 64 iterations of the GenASM-DC algorithm in parallel. Our hardware accelerator for GenASM-TB makes use of specialized per-PE TB-SRAMs, which eliminates the otherwise very high memory bandwidth consumption of traceback and enables e cient execution.
Technology-Level. With the help of 3D-stacked memory's vault-level parallelism, we can obtain 32× parallelism by performing 32 alignments in parallel in di erent vaults.
We have quantitatively evaluated three use cases of approximate string matching for genome sequence analysis (Section 10). We discuss four other potential use cases of GenASM, whose evaluation we leave for future work.
Read-to-Read Overlap Finding
Step of de Novo Assembly. De novo assembly [31] is an alternate genome sequencing approach that assembles an entire DNA sequence without the use of a reference genome. The rst step of de novo assembly is to nd read-to-read overlaps since the reference genome does not exist [152] . Pairwise read alignment (i.e., read-to-read alignment) is the last step of read-to-read overlap nding [102, 138] . As sequencing devices can introduce errors to the reads, read alignment in overlap nding also needs to take these errors into account. GenASM can be used for the pairwise read alignment step of overlap nding.
Hash- Table Based Indexing. In the indexing step of read mapping, the reference genome is indexed and stored as a hash table, whose keys are all possible xed-length substrings (i.e., seeds) and whose values are the locations of these seeds in the reference genome. This index structure is queried in the seeding step to nd the candidate matching locations of query reads. As we need to nd the locations of each seed in the reference text to form the index structure, GenASM can be used to generate the hash-table based index.
Whole Genome Alignment. Whole genome alignment [42, 136] is the method of aligning two genomes (from the same or di erent species) for predicting evolutionary or familial relationships between these genomes. In whole genome alignment, we need to align two very long sequences. Since GenASM can operate on arbitrary-length sequences as a result of our divide-and-conquer approach, whole genome alignment can be accelerated using the GenASM framework.
Generic Text Search. Although GenASM-DC is optimized for genomic sequences (i.e., DNA sequences), which are composed of only 4 characters (i.e., A, C, G and T), GenASM-DC can be extended to support larger alphabets, thus enabling generic text search. When generating the pattern bitmasks during the pre-processing step, the only change that is required is to generate bitmasks for the entire alphabet, instead of for only four characters. There is no change required to the edit distance calculation step.
As special cases of general text search, the alphabet can be de ned as RNA bases (i.e., A, C, G, U) for RNA sequences or as amino acids (i.e., A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V) for protein sequences. This enables GenASM to be used for RNA sequence alignment or protein sequence alignment [15, 16, 44, 67, 69, 90, 108, 126, 126, 128, 154, 157, 182] .
To our knowledge, this is the rst approximate string matching acceleration framework that enhances and acceler-ates the Bitap algorithm, and demonstrates the e ectiveness of the framework for multiple use cases in genome sequence analysis. Many previous works have attempted to improve (in software or in hardware) the performance of a single step of the genome sequence analysis pipeline. Recent acceleration works tend to follow one of two key directions [8] .
The rst approach is to build pre-alignment lters that use heuristics to rst check the di erences between two genomic sequences before using the computationally-expensive approximate string matching algorithms. Examples of such lters are the Adjacency Filter [177] that is implemented for standard CPUs, SHD [176] that uses SIMD-capable CPUs, and GRIM-Filter [91] that is built in 3D-stacked memory. Many works also exploit the large amounts of parallelism o ered by FPGA architectures for pre-alignment ltering, such as Gate-Keeper [10] , MAGNET [11] , Shouji [9] , and SneakySnake [13] . A recent work, GenCache [122] , proposes an in-cache accelerator to improve the ltering (i.e., seeding) mechanism of GenAx (for short reads) by using in-cache operations [1] and software modi cations.
The second approach is to use hardware accelerators for the computationally-expensive read alignment step. Examples of such hardware accelerators are RADAR [74] , FindeR [181] , and AligneR [180] , which make use of ReRAM based designs for faster FM-index search, or RAPID [65] and BioSEAL [88] , which target dynamic programming acceleration with processing-in-memory. Other read alignment acceleration works include SIMD-capable CPUs [38] , multicore CPUs [57, 109] , and specialized hardware accelerators such as GPUs (e.g., GSWABE [109] , CUDASW++ 3.0 [110] ), FPGAs (e.g., FPGASW [49] , ASAP [22] ), or ASICs (e.g., Darwin [162] and GenAx [55] ).
In contrast to GenASM, all of these prior works focus on accelerating only a single use case in genome sequence analysis, whereas GenASM is capable of accelerating at least three different use cases (i.e., read alignment, pre-alignment ltering, edit distance calculation) where approximate string matching is required.
We propose GenASM, an approximate string matching (ASM) acceleration framework for genome sequence analysis built upon our modi ed and enhanced Bitap algorithm. GenASM performs bitvector-based ASM, which can accelerate multiple steps of genome sequence analysis. We co-design our highly-parallel, scalable and memory-e cient algorithms with low-power and area-e cient hardware accelerators. We evaluate GenASM for three di erent use cases of ASM in genome sequence analysis for both short and long reads: read alignment, pre-alignment ltering, and edit distance calculation. We show that GenASM is signi cantly faster and more power-and area-e cient than state-of-the-art software and hardware tools for each of these use cases. We hope that GenASM inspires future work in co-designing algorithms and hardware together to create powerful frameworks that accelerate other bioinformatics workloads and emerging applications.