Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems
Vasimuddin Md Sanchit Misra Heng Li Srinivas Aluru
cs.DC cs.CE cs.PF q-bio.GN






Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems
    
Vasimuddin Md, Sanchit Misra, 
Parallel Computing Lab, 
 Intel Corporation, Bangalore, India. 
 Email: {vasimuddin.md,sanchit.misra}@intel.com 

Heng Li, 

Department of Biostatistics and Computational Biology, 
 Dana-Farber Cancer Institute, Boston, USA. 

Department of Biomedical Informatics, 
 Harvard Medical School, Boston, USA. 
 Email: hli@jimmy.harvard.edu  

Srinivas Aluru, 

School of Computational Science and Engineering, 
 Georgia Institute of Technology, Atlanta, USA. 
 Email: aluru@cc.gatech.edu 


    December 30, 2023
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================







Innovations in Next-Generation Sequencing are enabling generation of DNA sequence data at ever faster rates and at very low cost. For example, the Illumina NovaSeq 6000 sequencer can generate 6 Terabases of data in less than two days, sequencing nearly 20 Billion short DNA fragments called reads at the low cost of $1000 per human genome. Large sequencing centers typically employ hundreds of such systems. Such high-throughput and low-cost generation of data underscores the need for commensurate acceleration in downstream computational analysis of the sequencing data. A fundamental step in downstream analysis is mapping of the reads to a long reference DNA sequence, such as a reference human genome. Sequence mapping is a compute-intensive step that accounts for more than 30% of the overall time of the GATK (Genome Analysis ToolKit) best practices workflow. BWA-MEM is one of the most widely used tools for sequence mapping and has tens of thousands of users.

In this work, we focus on accelerating BWA-MEM through an efficient architecture aware implementation, while maintaining identical output. The volume of data requires distributed computing and is usually processed on clusters or cloud deployments with multicore processors usually being the platform of choice. Since the application can be easily parallelized across multiple sockets (even across distributed memory systems) by simply distributing the reads equally, we focus on  performance improvements on a single socket multicore processor. BWA-MEM run time is dominated by three kernels, collectively responsible for more than 85% of the overall compute time. We improved the performance of the three kernels by 1) using techniques to improve cache reuse, 2) simplifying the algorithms, 3) replacing many small memory allocations with a few large contiguous ones to improve hardware prefetching of data, 4) software prefetching of data, and 5) utilization of SIMD wherever applicable – and massive reorganization of the source code to enable these improvements. As a result, we achieved nearly 2×, 183×, and 8× speedups on the three kernels, respectively, resulting in up to 3.5× and 2.4× speedups on end-to-end compute time over the original BWA-MEM on single thread and single socket of Intel Xeon Skylake processor. To the best of our knowledge, this is the highest reported speedup over BWA-MEM (running on a single CPU) while using a single CPU or a single CPU-single GPGPU/FPGA combination.




Source-code: <https://github.com/bwa-mem2/bwa-mem2>





2019 IEEE.  Personal use of this material is permitted.  Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
 



§ INTRODUCTION





Innovations in Next-Generation Sequencing (NGS) have enabled generation of DNA sequence data at an enormous rate and very low cost. This is best exemplified by Illumina sequencers; a single Illumina NovaSeq 6000 can generate nearly 6 Terabases of data in a 44-hour run, sequencing nearly 20 billion short DNA fragments, called reads, of length 150 base-pairs each at the cost of only $1000 per human genome <cit.>. This has established short read sequencing using Illumina sequencers as the standard for modern genomics studies with large sequencing centers employing hundreds of such sequencers. 

The widespread use of NGS has enabled several applications – sequencing of individual genomes, sequencing genomes of new species, sequencing RNA samples to gather digital gene expression data, sequencing metagenomic samples of communities of microbial organisms, and single cell sequencing. Each of these usually requires mapping reads to one or more reference sequences or genomes. Sequence mapping is a compute intensive step; for example, it accounts for more than 30% of the overall time of the GATK (Genome Analysis ToolKit) best practices workflow <cit.>, one of the most prominent workflows for analyzing sequencing reads.



As the sequencers continue getting faster and cheaper at an exponential rate much faster than the Moore's law, we need commensurate speedups in sequence mapping. BWA-MEM <cit.> is one of the most popular tools used for mapping short reads to reference sequences as it is able to achieve both speed and accuracy. The software already has tens of thousands of users. The volume of data necessitates the use of distributed memory systems and most of it is processed on clusters and cloud deployments with multicore processors usually being the platform of choice.


In this work, we focus on accelerating BWA-MEM on multicore systems using efficient architecture-aware implementation. Since performance on multiple sockets can be achieved by just distributing the reads equally and load imbalance is usually not an issue, our efforts are focused on single socket performance. 

An important requirement is to maintain identical output. Many genomics studies can take years to complete. For these studies, it is critical that the output of the sequence mapping tool does not change over a long period of time. Thus, we maintain identical output while accelerating BWA-MEM to allow like-for-like replacement. This makes optimizing the tool extremely challenging. The compute time of BWA-MEM is dominated by three key kernels that collectively account for over 85% of the total compute time - (i) search for super maximal exact matches (SMEM) between reads and the reference sequence, (ii) suffix array lookup (SAL), and (iii) banded Smith Waterman (BSW) algorithm. There are several heuristics applied within the kernels and while connecting the kernels. The tool suffers from irregular memory accesses and irregular computation with too many branches and data dependencies making the use of SIMD extremely difficult. Attempts to make the computation or memory access more regular by modifying the algorithm usually result in change in output. 

Given the complexity of the BWA-MEM algorithm, only a few attempts have been made to accelerate it <cit.>. Most of these approaches use FPGAs or GPGPUs to accelerate one of the key kernels of BWA-MEM and report less than 2× speedup over the end-to-end compute time. One of them <cit.> demonstrates a 2.6× performance gain on end-to-end compute time by accelerating two out of the three kernels using four FPGAs.

In this paper, we present end-to-end optimization of BWA-MEM targeting multicore systems. We demonstrate the benefits of our improvements using an Skylake processor and an E5 v3 processor[Intel, Xeon and Intel Xeon Phi are trademarks of Intel Corporation or its subsidiaries in the U.S. and/or other countries. Other names and brands may be claimed as the property of others. Intel Corporation]
, but our optimizations are generic and can provide performance gains on any of the modern multicore processors.


Our key contributions are as follows.


    
  * We performed rigorous architecture-aware optimizations on the three kernels by 1) designing techniques to improve cache reuse, 2) simplifying the algorithms, 3) replacing many small fragmented memory allocations with a few large contiguous ones to improve hardware prefetching of data, 4) exploiting software prefetching of data, 5) utilizing SIMD wherever applicable, and 6) applying shared memory parallelism using OpenMP. We undertook a massive reorganization of the source code to enable these improvements.
    
  * We achieved 2×, 183× and 8× speedup, respectively, on SMEM, SAL and BSW kernels.
    
  * Our optimizations achieve up to 3.5× and 2.4× speedups on end-to-end compute time over the original BWA-MEM on single thread and single socket, respectively, of Skylake processor.
    
  * Since our output is identical to the original BWA-MEM, it can seamlessly replace the original.
    
  * We present detailed performance analysis of our implementation, and use the resulting insights to conduct architecture bottleneck analysis when using the Skylake processor for BWA-MEM.


To the best of our knowledge, this is the highest reported speedup over BWA-MEM while using a single CPU or a single CPU-single GPGPU/FPGA combination. 








§ A BRIEF OVERVIEW OF BWA-MEM



In what follows, we use capital letters X, Y, etc. to denote DNA sequences, modeled as strings over the alphabet Σ = {A,C,G,T} representing the four possible nucleotides (also called bases). The bases within a sequence are denoted using small letters like a, b, c, etc. We use |X| to represent the length of sequence X. Let X[i] denote the base at position i, and X[i,j] (i≤ j) denote the substring X[i]X[i+1]X[i+2]⋯ X[j]. The bases `A' and `T' are considered complements of each other and `C' and `G' are complements of each other. We denote the complement of a base c as c̅.



 §.§ Short Read Sequence Mapping Problem


Let Q be a short read and R be a long reference sequence. The sequence mapping problem is to find the best matches of Q in R. For short reads, typical lengths of Q range from 50 to 250. |R| can range anywhere from a few hundred thousand (for bacterial genome) to several billions (for plant genomes). An important special case is the human genome, at 3 billion base-pairs long.

Most popular sequence mapping tools, including BWA-MEM, use the seed and extend method.
In the seeding phase, matches for short substrings of the reads (called seeds) are found in the reference sequence while allowing for zero or very small number of differences in the match. This gives the potential locations where the entire read can match. In the extension phase, the seed matches are extended on both sides to cover the entire query, typically using a dynamic programming based algorithm, and the best matches are reported. BWA-MEM uses an FM-index of the reference sequence to perform seeding.







 §.§ FM-index





The FM-index of a reference sequence is based on its Burrows-Wheeler Transform (BWT) and suffix array. Figure <ref> shows how to construct the BWT, suffix array and FM-index of an example sequence R. First, R is appended with the character $∉Σ which is lexicographically smaller than all characters in Σ. Subsequently, all the rotations of R are obtained (Rotations(R)). Lexicographically sorting the rotations gives the BW matrix. The BWT (B) is the last column of this matrix.
The suffix array (S) stores the original positions in R of the first bases of these rotations, equivalent to the sorted order of the suffixes of R.
All the exact matches of a query are prefixes of the rotations in the BW-matrix. These matches are located in contiguous rows of the BW-matrix since the rotations are lexicographically sorted. Therefore, all the matches of a query can be represented as a range of rows of the BW-matrix. This range is called the SA interval of the query. For example, in Figure <ref>, matches of sequence “AC” correspond to the SA interval [1,2].

To facilitate fast search of BWT, FM-index is used  <cit.>. It consists of D and O data structures. D[x] is the number of bases in R[0,|R|-1] (excluding $) that are lexicographically smaller than x ∈Σ, and O[x,i] is the number of occurrences of x in B[0,i].



 §.§ BWA-MEM algorithm



The BWA-MEM algorithm consists of the following steps.


    
   SMEM: Perform seeding by searching for super maximal exact matches (SMEM) between the read and the reference sequence. Formally, maximal exact matches (MEMs) are exact matches between substrings of two sequences that cannot be further extended in either direction. An SMEM is a MEM that is not contained in any other MEMs on the query sequence <cit.>. BWA-MEM uses FM-index to search for SMEMs and outputs SA intervals of the SMEMs.
    
   SAL: Suffix array lookup is performed to get the coordinates in the reference sequence from the SA intervals.
    
   CHAIN: BWA-MEM chains the seeds that are collinear and close to each other. The chains thus formed are used to filter seeds for alignment.
    
   BSW: The seeds are extended using dynamic programming (DP) based banded Smith-Waterman (BSW) alignment algorithm that computes only a diagonal band of the DP matrix. The diagonal band may shrink as the computation moves from the top row to the bottom row.
    
    
    
    
    

    
   SAM-FORM: The algorithm concludes by formatting the alignment output in the SAM format.
    





 §.§ Break up of computation time spent



























To identify compute hotspots for acceleration,
we performed rigorous run-time profiling of BWA-MEM (version 0.7.25) by hand instrumenting the source code with trackers. Table <ref> reports percentage of the run-time spent in various steps of BWA-MEM. 



Clearly, SMEM, SAL and BSW are the most time consuming steps accounting for 86.5% and 85.7% of the total run-time of BWA-MEM, on D1 and D4 datasets, respectively. Therefore, we select them as targets for acceleration. 




 §.§ Potential for improvement of the key kernels








In this section, we analyze the performance characteristics of the three kernels of BWA-MEM to identify target areas for performance improvement. To do so, we extracted the source code corresponding to each kernel from BWA-MEM and created a benchmark for each kernel. To prepare input datasets for the benchmarks, we executed BWA-MEM using read datasets and intercepted inputs to each of the kernels. The benchmarks produce the exact same output as the corresponding kernels in the original BWA-MEM.



  §.§.§ SMEM


As shown in Section  <ref>, Table  <ref>, the SMEM kernel executes nearly 17 Billion instructions for 60,000 reads; that is; nearly 285 thousand instructions per read. The main reason for this is that BWA-MEM uses a compressed version of FM-index. Due to excessive compression, it requires a large number of instructions to compute the SA interval. Moreover, the kernel suffers from high LLC miss rate resulting in a significantly high average memory latency (24).



  §.§.§ SAL


SAL kernel just retrieves the coordinate of reference sequence from the SA interval. What should have been a simple array lookup requires nearly 5000 instructions per lookup due to the use of compressed suffix array.



  §.§.§ BSW


Due to the irregular nature of the BSW algorithm with significant branching and short loops, it uses a scalar implementation. As we explain in Section <ref>, BSW is a compute bound kernel. So, while memory access is not an issue, BSW implementation of BWA-MEM is instruction bound, due to the scalar implementation, resulting in lower performance.




§ MODIFICATIONS APPLIED TO ENTIRE CODE





 §.§ Reorganization of the Workflow






Original BWA-MEM processes a chunk of reads at a time. Figure <ref> shows how we transform the workflow of processing one chunk. For each read in a chunk, the original implementation processes it through all the steps until BSW before moving on to the next read. It maintains the BSW output for the reads until the entire chunk is processed. The BSW output of each read is then converted to SAM format. Multithreading is done using pthreads by dynamically distributing the reads across threads. On the other hand, we divide the reads in a chunk into batches. Each step is executed on all the reads in a batch before moving to the next step. Multithreading is done using OpenMP by dynamically distributing the batches across threads.

The benefit of such re-organization is that it enables the use of SIMD parallelism across different reads in a batch. We make use of this for BSW.



 §.§ Improving Inefficient Memory Allocation


The BWA-MEM implementation allocates and de-allocates small blocks of memory frequently, which is expensive, does not provide hardware prefetcher with a long predictable pattern of memory access, and does not allow cache reuse. Instead, we allocate all the required memory once in large contiguous blocks and reuse it across the batches. Contiguous memory allocation provides benefits of hardware prefetching, reducing the memory latency. Moreover, since same buffer is used repeatedly across batches, it also improves cache reuse.




§ OPTIMIZATIONS APPLIED TO SMEM AND SAL





 §.§ FM-Index Structure



pseudocode


The implementation of FM-index in BWA-MEM differs
from the description in Section <ref> in two ways. First, SMEM algorithm requires the FM-index of the sequence formed by concatenating the reference sequence with its reverse complement. Second, in order to reduce the memory footprint, the FM-index is compressed by dividing O into buckets of size η. The compressed data structure O_c has only |R|/η entries. In each entry i, two items are stored – 1) count: the count values of the bases up to bucket i; thus, for each base c, O_c[i].count(c) holds the counts O[c, i ×η], and 2) bwt: the substring of BWT, B[i ×η,(i + 1) ×η - 1]. The O values of any position can be computed using the O_c count and bwt (Algorithm  <ref>). A power-of-2 number is picked for the value of η to replace the expensive division and modulo operations (lines 1 and 2) with right-shift and bitwise AND operations. Base c and bwt use 2-bit representation using { 0,1,2,3 } for {A,C,G,T}. 

The count of base c in the first y positions of bwt (line 4) is obtained by using bitwise operations.





 §.§ SMEM Search Algorithm



pseudocode


pseudocode


pseudocode


Algorithms  <ref>,  <ref> and  <ref> provide the details of SMEM search algorithm as given in  <cit.>. The dots signify details of the algorithm that are omitted as they are not relevant to this discussion. The software prefetching, introduced by us, is discussed in the next Section. Algorithm  <ref> returns all the SMEMs passing through position i_0 in X. It starts at position i_0 and finds all the extensions in the forward direction that have matches. For all matches, it maintains the bi-interval, (k,l,s), where - (i) k is the starting position of the SA interval of the matching substring of X, (ii) l is the starting position of the SA interval of the reverse complement of the matching substring, and (iii) s is the size of the SA interval. For all the forward extensions, it finds the longest backward extensions that have matches to find all the SMEMs. The forward extension algorithm calls backward extension, which in turn calls Get_O that accesses O_c. The successive inputs (k,l,s) to backward extension have no particular order or relationship. Therefore. access to O_c has no spacial or temporal locality and is completely irregular. Given the huge size of O_c, this leads to frequent cache misses making SMEM kernel memory latency bound. 

While the compression of FM-index reduces memory footprint, there is also an additional benefit. More accesses fall into the same O_c entry, improving data locality. In particular, as the matching string gets longer, the corresponding SA interval gets shorter and the likelihood of both k and k+s entries falling in the same cache line increases. With compressed FM-index, this should happen more frequently. The trade off is the increase in runtime for processing a base from O(1) to O(h).




 §.§ Applying Software Prefetching to SMEM Kernel



Given the irregular memory access to O_c, hardware prefetching will be ineffective. However, note that each access to O_c results in new values of (k,l,s) that, if a match, decides the memory locations of a later access to O_c. We can prefetch this location using software prefetching. Therefore, whenever new values of (k,l,s) are computed that are likely to be used in the future to access O_c, we software prefetch the corresponding locations of O_c. During forward extension, if a new match ((k,l,s)) is found, it is either used in the next forward extension or appended to Curr which is used for a future backward extension. Therefore, we prefetch the corresponding locations of O_c. Similarly, during backward extension, any new match ((k,l,s)) that is added to Curr is used for a future backward extension and, thus, we prefetch the corresponding locations of O_c.
However, software prefetching will only be useful if we can hide the latency of prefetching by computation. When a match found during forward extension is used for the next forward extension, or when the number of matches in Prev are small during backward extension, there may not be sufficient compute to hide memory latency. Since backward extension stage dominates the run time and Prev is sufficiently large frequently enough, software prefetching is quite useful. But it can not alleviate memory latency completely.

In order to get sufficiently large compute to hide the memory latency for every O_c access, we also tried processing multiple queries simultaneously in a round robin fashion performing one call to backward extension (or forward extension) per query every time and software prefetching according to the new match produced. However, given the intricacy of the algorithm and number of different branches possible, this approach resulted in a much larger number of instructions. The additional advantage of alleviating memory latency completely could not overcome the disadvantage of extra instructions.



 §.§ Choosing Compression Factor


For our implementation, we choose the value of η to be 32 due to the following reasons.
It is better to ensure that the size of one entry of O_c is less than that of one cache line. Otherwise, we will need to prefetch multiple cache lines per O_c access that can cause more data to be read from memory resulting in higher memory bandwidth requirement. 

Each count needs 4 bytes (unsigned integer). Thus, in a O_c entry, 16 bytes are consumed to store the count values. That leaves a maximum of 48 bytes for the BWT substring. We use power-of-2 value for η to prevent expensive division and modulo operations. Thus, we can only use 32 bytes for the BWT substring. BWA-MEM uses the value of η as 128 and uses 2-bits to represent each base in the BWT string, thus consuming 32 bytes for the BWT string. However, extracting the count of a particular base from a BWT string of length 128 requires a large number of instructions. Moreover, 2-bit representation requires significant bit manipulation making it expensive. Therefore, we opt for a one byte representation of bases to enable partial vectorization of the  computation of occurrences, thus choosing the value of η as 32. We perform a byte level compare using AVX2 to get a 32-bit mask containing 1 for match and 0 for mismatch. Consequently, we use a 32-bit popcnt instruction on the mask to get the count. We store each entry of O_c in a cache aligned location for effective prefetching by adding padding of 16 bytes to use a full cache line per entry.



 §.§ Improvements to Suffix Array Lookup



Given a position i in the suffix array, S, SAL returns the corresponding coordinate in reference sequence, j, as shown in Equation <ref>.


    j = S[i]


Similar to FM-index, BWA-MEM stores a compressed version of S with compression factor of 128. It obtains the value of S[i] through an intricate algorithm requiring access to multiple locations in the compressed S and FM-index. As a result, it requires nearly 5000 instructions on an average for each S[i] and a large number of memory accesses.

There is sufficient memory available to store the suffix array for most reference sequence sizes. For example, for entire human genome, we need about 48 GB of memory. Therefore, we use the uncompressed SA and use the expression in Equation <ref> to get the coordinate of reference sequence.




§ OPTIMIZATIONS APPLIED TO BSW






















 §.§ Banded Smith Waterman Algorithm Used in BWA-MEM



The BSW algorithm used in BWA-MEM is derived from the Smith Waterman (SW) algorithm, which is a dynamic programming (DP) based sequence alignment algorithm. Given a pair of DNA sequences X and Y, both SW and BSW compute the cells of a DP matrix H of size |X| × |Y| in a row-wise manner. Here, the value H[i,j] at each cell represents the alignment score between sub-sequences X[0,i-1] and Y[0,j-1]. The key difference in the computation pattern is that, while SW computes all the cells of the DP matrix, 
[label=(*)]
    
  * BSW restricts the cell computations to within a certain band size (λ) around the main diagonal, 
    
  * BSW aborts the matrix computations, if all the scores in a row are zero or the best score of the current row drops by a certain threshold from best score so far, and
    
  * after computing a row, BSW adjusts the value of λ based on number of cells with zero value from both ends.

The top matrix of Figure <ref> illustrates the full matrix and a diagonal band.


In BSW, computation of a cell (i,j) of H uses the following recursion.



    H[i,j] = max{ H[i-1,j-1] + S(X[i],Y[j]), E[i,j], F[i,j] }
       E[i+1,j] = max{ H[i,j] - 𝗀_o, E[i,j] - 𝗀_e, 0 }
       F[i,j+1] = max{ H[i,j] - 𝗀_o, F[i,j] - 𝗀_e, 0 }


Where, 𝗀_o and 𝗀_e are gap open and extension penalties.
S(a,b) = m (match score), if a=b, else S(a,b) = m' (mismatch score). For 1 ≤ i ≤ |X|, 1 ≤ j ≤ |Y|, the matrices are initialized as,







 














    H[i,0]    =  max{ 0, H[0,0] - (𝗀_o + 𝗀_e . i) }, 
    
        H[0,j]    = max{ 0, H[0,0] - (𝗀_o + 𝗀_e . j) },  
    
        E[i,1]    = F[1,j] = 0


Since BSW extends a seed, the initial matrix score H[0,0] is the match score of the seed. Given the computation of a cell is only dependent on the three neighbors, BSW in BWA-MEM only maintains one row each of E, F, and H. Smaller sequence sizes ensure that these arrays are present in caches and only the sequence pairs need to be fetched from memory, making BSW a compute bound problem.



 §.§ Vectorization approaches and challenges


Given that BSW is compute bound, we chose to vectorize it using SIMD parallelism. Vectorization of standard SW has been studied extensively. There are two popular approaches: intra-task  and inter-task. We see the following challenges in application of these approaches to BSW.



  §.§.§ Intra-task


Intra-task approach vectorizes the computation of SW of one pair. It either leverages the independence between anti-diagonal cells of the matrix to vectorize across them or uses striped SIMD approach. However, this kind of approach is unwieldy for BSW. Short sequences, computation only within the band, and potential reduction in the band size due to band adjustment, restrict the available parallelism for vectorization.



  §.§.§ Inter-task


In inter-task approach, each vector lane computes the matrix for different sequence pairs. The irregularity of matrix sizes, number of cells computed, and the positions within matrices of cells computed across matrices, makes extracting performance using inter-task approach difficult. 

Moreover, in BWA-MEM, whether a seed is extended using BSW or not depends on the previously extended seeds of a read. Such dependency restricts the parallelism across seeds of a read. 





Due to availability of better parallelism, we utilize inter-task vectorization for accelerating BSW (Figure <ref>). 





 §.§ Inter-Task Vectorization



Our approach for inter-task vectorization is illustrated in Figure <ref> and is based on our previous work <cit.>. Let W be the SIMD width. We process W sequence pairs at a time. For each cell (i,j) that is computed, each vector lane computes the corresponding cell (i,j) for all matrices. Therefore, even if one of the sequence pairs requires a particular matrix cell to be computed, the corresponding matrix cell is computed for all the sequence pairs leading to wasteful cell computations.



  §.§.§ Sorting

While any aborted sequence pair can be replaced with a new one to avoid idle vector lanes, the overhead of the operation supersedes its benefits.
Wasteful cell computations due to variations in sequence lengths can be better curbed by ensuring uniformity in sequence lengths. We use radix sort to sort the tasks by their respective sequence lengths, and then group together tasks with the same or close sequence lengths to ensure uniformity of tasks filling vector lanes.



  §.§.§ Increasing parallelism


Given the dependency between seeds, the parallelism is only available across reads. However, the count of seeds per read can vary a lot leading to imbalance of computation across reads. A dynamic allocation of reads to vector lanes can potentially balance computation across vector lanes. However, given the imbalance in number of computed seeds across reads, it requires a large number of reads in a batch to succeed. Large batches are not possible due to memory constraints as we have to maintain the metadata corresponding to each read between steps of the BWA-MEM algorithm. 
Therefore, we manage the dependency between the seeds by first extending all the seeds of a read, and then post process them to filter out the ones that should not have been extended.
 


  §.§.§ AoS to SoA conversion of sequence pairs


We convert the input sequences from AoS to SoA format to allow reading of corresponding bases from W sequence pairs using vector load operation instead of a gather operation.



 §.§ Overview of operations


For a given row, we do the following computations. 
(a) To perform cell computations within a given band, we utilize instructions such as cmp, blend, max, mov, add, and sub. Similar instructions are used to identify the best score and its position.
(b) To find the global score, at each cell in a row, we check whether the end of query is reached for all the pairs using cmp instruction and update the global score and its position using cmp and blend instructions.

(c) After the row computations, the cell range (for next row) is adjusted again based on the number of empty cells from both ends of the row.
From both ends of the row, with one cell at a time, we use cmp and blend vector instructions to check for empty cells and update the range.
(d) Moreover, we use mask and cmp instructions for maintaining the correct values for aborted sequence pairs.



  §.§.§ Precision


Smaller integer-widths provide more vector lanes.
Since matrix scores are proportional to sequence lengths, it is possible to use different integer-widths for matrix computation depending on sequence lengths. In BSW optimization, we use 8-bit or 16-bit implementations depending on the sequence length.








§ RESULTS



Software and workloads used in performance tests may have been optimized for performance only on Intel microprocessors. Performance tests, such as SYSmark and MobileMark, are measured using specific computer systems, components, software, operations and functions. Any change to any of those factors may cause the results to vary. You should consult other information and performance tests to assist you in fully evaluating your contemplated purchases, including the performance of that product when combined with other products. For more information go to www.intel.com/benchmarks.

Benchmark results were obtained prior to implementation of recent software patches and firmware updates intended to address exploits referred to as "Spectre" and "Meltdown".  Implementation of these updates may make these results inapplicable to your device or system.




 §.§ Experimental Setup




  §.§.§ System Configuration






We evaluated our solution on the latest Intel Xeon Scalable family of processor (Skylake). Our specific Skylake processor is detailed in Table  <ref> and referred to as SKX from here on. Each AVX512 VPU can process multiple 8-bit and 16-bit integers thus capable of SIMD widths of 64 and 32, respectively. We used 2018 for performance analysis by measuring hardware performance counters.


While our implementation achieves best performance by running on latest processors with AVX512 support from heavy use of AVX512 instructions, we also support alternate AVX2-only and scalar-only executions for older generation multicore processors. To confirm the wide applicability of our improvements, we evaluated our implementations on two generations of processors – the latest Intel Xeon Scalable family of processors (Skylake) and Intel E5 v3 family (Haswell). Our specific Skylake and Haswell processors are detailed in Table  <ref> and referred to as SKX and HSW, respectively, from here on. Each AVX512 VPU can process multiple 8-bit and 16-bit integers thus capable of SIMD widths of 64 and 32, respectively; similarly, each AVX2 VPU can process multiple 8-bit and 16-bit integers thus capable of SIMD widths of 32 and 16, respectively. We used 2018 for performance analysis by measuring hardware performance counters.

We perform all our experiments on a single sockets of SKX (28 cores) and HSW (18 cores) and use numactl utility to force all memory allocations to one socket. For multi-threaded runs, we run 2 threads per core to get the benefit of hyper threads. Since optimizing file IO is beyond the scope of this paper, we do not include file IO time in any of our results.



  §.§.§ Datasets

To demonstrate the performance of our improvements, we make use of the following real datasets. We use the first half of the human genome (version Hg38), containing nearly 1.5 billion nucleotides, as the reference sequence for all our experiments. We use five different read datasets covering all prominent short read lengths (Table <ref>) to demonstrate the versatility of our implementation. We received the read datasets D1 and D2 directly from the Broad Institute to evaluate the performance of our implementation.





  §.§.§ Correctness


For each experiment, we verified that our output is exactly identical to that of original BWA-MEM.




 §.§ Performance Evaluation of Key Kernels On a Single Thread

In this section, we evaluate the benefits of our improvements on the individual key kernels using real datasets on SKX. We extracted the kernels from both the implementations so that we are comparing only the kernels and nothing else.




  §.§.§ SMEM search using FM-index











Table  <ref> evaluates our improvements to the SMEM kernel. Our modifications in the design to use the value of bucket size η as 32 and vectorization to reduce the number of instructions has resulted in a greater than 2× reduction in the total number of instructions and the number of load instructions. 

However, lower value of η results in less values of k and k+s falling into the same bucket, thus resulting in lower cache reuse. This is shown in the increase in the number of LLC misses and average latency of memory accesses between “Original” and “Optimized minus S/W prefetching” implementations. The use of software prefetching dramatically reduces the number of LLC misses (by nearly 3×) and average latency of memory accesses at the cost of slightly higher instruction count. The change in performance counters is clearly reflected in the time consumed. Overall, our improvements achieved 2× performance gain over the original.




  §.§.§ Suffix Array Lookup (SAL)







Table  <ref> shows the benefit of simplifying the SAL kernel. There is a dramatic reduction in the number of instructions per SA offset (nearly 200×). Therefore, despite the increase in number of LLC misses per load, and thereby average memory latency, we see a speedup of nearly 183×. There is further scope of performance improvement by using software pre-fetching. However, for our optimized implementation the time consumed in this kernel is negligible compared to the full application run time. Therefore, there isn't much benefit in pursuing it.




  §.§.§ Banded Smith-Waterman



















Table <ref> shows the performance benefit of our vectorized implementations at 8-bit and 16-bit precisions with or without sorting. For this experiment, we only used the sequence pairs for which 8-bit precision was sufficient. The results show clear benefit of sorting which provides 1.5× to 1.7×  performance boost by increasing the chances of pairs that are processed simultaneously using SIMD parallelism having similar lengths. Overall, our 8-bit and 16-bit implementations achieve nearly 11.6× and 6.7× speedup, respectively, over the original implementation.




Table  <ref> uses hardware performance counters to explain the speedup. BSW is a compute-intensive task.

Consequently, any reduction in number of instructions executed would result in reduction in execution time. 

Vectorization led to reduction in the number of instructions executed by 13.85× (for 8-bit implementation) over the original implementation of the BSW kernel that was a scalar implementation. The reduction in IPC is explained by the fact that a large majority of the instructions in our optimized code are SIMD instructions. Therefore, with two ports for SIMD instructions, an IPC of 2.17 is expected while accounting for a few scalar instructions. In contrast, there are 4 ports for scalar ALU instructions explaining an IPC of 3.17 for the original implementation. 




Note that 16-bit and 8-bit integer precision values provide SIMD widths of 32 and 64, respectively, with AVX512 facilitating potential speedups of 32× and 64×, respectively. In order to understand the gap between these ideal speedups and achieved speedups, Table  <ref> presents the breakdown of time spent in different components of the optimized 8-bit BSW. 
Run-time profiling shows that only 43% of the total time is spent in calculating the DP matrix; while remaining time is spent in pre-processing of the input sequences (which converts the input sequences from AoS to SoA order), and cell range adjustment. AoS to SoA conversion is performed specifically for the vectorized implementation and is not needed for the original scalar implementation.


Owing to multiple factors such as variations in sequence lengths and aborting of matrix computations (discussed in section <ref>), our inter-task vectorization ends up computing wasteful cells.
Empirical analysis shows that the useful cells are roughly half of the total cells computed. Therefore, only 21.5% of the run time is spent in computation of useful cells of the DP matrices.








 §.§ End-to-end Performance Evaluation on Single Socket










In this section, we present the comparison of the end-to-end compute time of the original and our optimized implementation of BWA-MEM, with all the optimized kernels integrated.



  §.§.§ Scaling with respect to the number of cores








We evaluate the multicore scaling of our implementation and three integrated kernels in (Figure <ref>). All the kernels demonstrate good scaling from 1 core to 28 cores on single socket closely matching the scaling of the original BWA-MEM with our optimized implementation achieving better scaling for the memory intensive SMEM and SAL kernels for the D5 dataset. We see a dip in scaling performance of the overall application for our implementation (D1: 22×, D5: 20×). The decline in scaling stems from the part of the application other than the three kernels. That is the part of the implementation that has not been optimized yet. Though all three key kernels achieve greater than 25× scaling, some of the components in Misc are memory bandwidth bound and do not scale beyond 15×.




  §.§.§ Time to solution












We compared the time to solution of our optimized implementation with the original implementation using five different real datasets (Table <ref>). 

The run-time of the two implementations, and speedup of our implementation over original BWA-MEM, are depicted in Figure <ref> (single thread) and Figure <ref> (single socket). 


The run-time of our two implementations and original BWA-MEM, and speedup of our implementations over original BWA-MEM, are depicted in Figures <ref> & <ref> (single thread) and Figures <ref> & <ref> (single socket). 



All the three kernels sustain their benchmark speedup in the application. In fact, SAL kernel's run-time barely contributes to the overall run time due to the 183× speedup. However, application performance benefit for the BSW kernel is lowered owing to the processing of extra seeds (Section  <ref>). On an average, we see 14% extra sequence pairs aligned by our method. 
For dataset D2, we observed that there were 13.5% extra seed pairs, but resulted in nearly 1.43× more BSW time, incurring a loss of speedup over BSW of BWA-MEM by 1.47 ×.


On overall application, we see speedup range of 2.6× - 3.5× on single thread, and of 1.7× - 2.4× on single socket.


For dataset D2 ,on SKX, we observed that there were 13.5% extra seed pairs, but resulted in nearly 1.43× more BSW time, incurring a loss of speedup over BSW of BWA-MEM by 1.47 ×.


On overall application: on SKX, we see speedup range of 2.6× - 3.5× on single thread, and of 1.7× - 2.4× on single socket; while, we see speedup range of 2.3× - 3.0× on single thread, and of 1.9× - 2.7× on single socket, on HSW.














































§ RELATED WORK




Owing to the complexity of BWA-MEM, a majority of the attempts to speed it up are confined to accelerating one of the three kernels, either on GPGPU or on FPGA <cit.>. 

Additionally, these methods pipeline the non-optimized kernels on the host CPU for execution. Out of the three kernels, given that standard SW has been extensively targeted for acceleration and BSW of BWA-MEM is a variant of standard SW, a lot of studies <cit.> have targeted only BSW for acceleration. Due to dependencies between the seeds within a read, all of the mentioned acceleration approaches for BSW kernel employ intra-task parallelism.
Excluding  <cit.>, the other approaches manage to achieve 1.6× to 3× speedup on BSW and 1.45× to 2× overall. 
Chang et al. <cit.> exclusively target SMEM  optimizations on FPGAs and demonstrate 4× and 1.26× performance improvement on SMEM and overall application, respectively. Due to  pipelining of computations between host CPU and accelerator, the benefits of BSW speedup are contingent on speed of non-optimized kernels. 
Ahmed et al. <cit.> target all the three kernels for optimizations. They use 4 FPGAs to speedup SAL and BSW kernels by 2.8× and 5.7× respectively. Also, they apply algorithmic enhancements to optimize SMEM kernel on host CPU, speeding it up by 1.7× and overall performance by 2.6×.
However, their output may differ from BWA-MEM output, due to bypassing of some of the heuristics used in BWA-MEM BSW. To the best of our knowledge, no published work contains a holistic architecture-aware optimization of BWA-MEM software on multicore systems.








FM-index methods and SW based alignment are widely used techniques in NGS data analysis and have been classified as computational building blocks <cit.>.
Due to the intricacy of FM-index based algorithms, few works considered optimization of these problems <cit.> and most of them have targeted the simpler, exact search algorithm (find exact matches of full query sequences).
Chacon et al. <cit.> developed a clever technique for exact search, which reduces the number of memory accesses. Application of their technique to SMEM search is non-trivial. Additionally, they tried software prefetching for exact search <cit.> with minimal benefit due to an inefficient design. 

Misra et al. <cit.> demonstrated significant performance gains on exact search with software prefetching by batching multiple queries together to hide memory latency.




Intra- and inter-task vectorization are both popular approaches for SW acceleration <cit.>. However, BSW of BWA-MEM has some key variations from standard SW as discussed in section <ref>. Therefore, there have been only a few attempts to accelerate exact BSW kernel from BWA-MEM.





§ CONCLUSION AND FUTURE WORK




In this paper, we presented an efficient implementation of BWA-MEM for multicore systems while maintaining identical output so that the current users of BWA-MEM can seamlessly switch to the new implementation. We presented our improvements through architecture aware optimization to speedup the three key kernels - SMEM, SAL and BSW - by 2×, 183×, and 8×, respectively. This results in an overall performance gain of 3.5× and 2.4× over the compute time of BWA-MEM on a single thread and single socket, respectively, of an Intel Skylake processor. 


Our implementation has been open sourced so that the users of BWA-MEM can benefit from increased performance.

We were successful in making the SAL kernel virtually insignificant contributor to the overall run time. The other two kernels are instruction bound and SMEM is also partly memory latency bound. Given the irregular structure of the algorithms, SMEM can not benefit from vectorization, while gains for BSW due to vectorization are limited. Therefore, processors that rely heavily on SIMD performance can only achieve limited performance on BWA-MEM. While better support for gather operations will help, an architecture that does not rely on SIMD for performance will be better.

As future work, we plan to make the following improvements - (i) 
attempt to improve the parts of BWA-MEM apart from the three kernels that we focused on in this paper in order to improve multicore scaling and performance of the overall application further, and (iii) explore if there is additional scope of reducing average memory latency of SMEM kernel and the instruction counts of SMEM and BSW kernels. 


IEEEtran



