ApHMM: Accelerating Profile Hidden Markov Models for Fast and Energy-efficient Genome Analysis

2.1 Profile Hidden Markov Models (pHMMs)

High-level Overview. Figure 1 shows the traditional structure of pHMMs. pHMMs represent a sequence or a group of sequences using a certain graph structure with a fixed number of nodes for every character of represented sequences. Visiting nodes, called states, via directed edges, called transitions, are associated with probabilities to identify differences at any position between the represented sequences and other sequences. States emit one of the characters from the defined alphabet of the biological sequence (e.g., A, C, T, and G in DNA sequences) with a certain probability. Transitions preserve the correct order of the represented sequences and allow making modifications to these sequences. We explain the detailed structure of pHMMs in Supplemental Section 1.

2.2 Baum-Welch Algorithm

The probabilities associated with transitions and states are essential for identifying similarities between sequences. The Baum-Welch algorithm provides a set of equations to update and use these probabilities accurately. To calculate the similarity score of input observation sequences in a pHMM graph, the Baum-Welch algorithm [89] solves an expectation-maximization problem [98, 99, 100, 101], where the expectation step calculates the statistical values based on an input sequence to train the probabilities of pHMMs. To this end, the algorithm performs the expectation-maximization based on an observation sequence S for the pHMM graph \(G(V,A)\) in three steps: (1) forward calculation, (2) backward calculation, and (3) parameter updates.

Forward Calculation. The goal of the forward calculation is to compute the probability of observing sequence S when we compare it with the sequence \(S_{G}\) that the pHMM graph \(G(V,A)\) represents. Equation (1) shows the calculation of the forward value \(F_{t}(i)\) of state \(v_i\) for character \(S[t]\). The forward value, \(F_{t}(i)\), represents the likelihood of emitting the character \(S[t]\) at position t of S in state \(v_{i}\) given that all previous characters \(S[1 \dots t-1]\) are emitted by following an unknown path forward that leads to state \(v_{i}\). \(F_{t}(i)\) is calculated for all states \(v_i \in V\) and for all characters of S. Although t represents the position of the character of S, we use the timestamp term for t for the remainder of this article. To represent transition and emission probabilities, we use the \(\alpha _{ji}\) and \(e_{S[t]}(v_{i})\) notations as we define in Supplemental Section ??: (1) \(\begin{equation} F_{t}(i) = \sum _{j \in V} F_{t-1}(j) \alpha _{ji} e_{S[t]}(v_{i}) \hspace{5.0pt}i \in V, \hspace{5.0pt}1 \lt t \le n_{S}. \end{equation}\)

Backward Calculation. The goal of the backward calculation is to compute the probability of observing sequence S when we compare S and \(S_G\) from their last characters to the first characters. Equation (2) shows the calculation of the backward value \(B_{t}(i)\) of state \(v_i\) for character \(S[t]\). The backward value, \(B_{t}(i)\), represents the likelihood of emitting \(S[t]\) in state \(v_{i}\) given that all further characters \(S[t+1 \dots n_{S}]\) are emitted by following an unknown path backwards (i.e., taking transitions in reverse order). \(B_{t}(i)\) is calculated for all states \(v_i \in V\) and for all characters of S: (2) \(\begin{equation} B_{t}(i) = \sum _{j \in V} B_{t+1}(j) \alpha _{ij} e_{S[t+1]}(v_{j}) \hspace{5.0pt}i \in V, ~1 \le t \lt n_{S}. \end{equation}\)

Parameter Updates. The Baum-Welch algorithm uses the values that the forward and backward calculations generate for the observation sequence S to update the emission and transition probabilities in \(G(V, A)\). The parameter update procedure maximizes the similarity score of S in \(G(V, A)\). This procedure updates the parameters shown in Equations (3) and (4). The special \([S[t] = X]\) notation in Equation (4) is a conditional variable such that the variable returns 1 if the character X matches with the character \(S[t]\) and returns 0 otherwise: (3) \(\begin{equation} \alpha ^{*}_{ij} = \dfrac{\sum \nolimits _{t=1}^{n_{S}-1} \alpha _{ij}e_{S[t+1]}(v_{j})F_{t}(i)B_{t+1}(j)}{\sum \nolimits _{t=1}^{n_{S}-1}\sum \nolimits _{x \in V} \alpha _{ix}e_{S[t+1]}(v_{x})F_{t}(i)B_{t+1}(x)} \quad \forall \alpha _{ij} \in A, \end{equation}\) (4) \(\begin{equation} e^{*}_{X}(v_{i}) = \dfrac{\sum \nolimits _{t=1}^{n_{S}} F_{t}(i)B_t(i)[S[t] = X]}{\sum \nolimits _{t=1}^{n_{S}} F_{t}(i)B_t(i)} \quad \forall X \in \Sigma , \forall i \in V. \end{equation}\)

2.3 Use Cases of Profile HMMs

Error Correction. The goal of error correction is to locate the erroneous parts in DNA or genome sequences and replace these parts with more reliable sequences [102, 103, 104, 105, 106, 107] to enable more accurate genome analysis (e.g., read mapping and genome assembly). Apollo [24] is a recent error correction algorithm that takes an assembly sequence and a set of reads as input to correct the errors in an assembly. Apollo constructs a pHMM graph for an assembly sequence to correct the errors in two steps: (1) training and (2) inference. First, to correct erroneous parts in an assembly, Apollo uses reads as observations to train the pHMM graph with the Baum-Welch algorithm. Second, Apollo uses the Viterbi algorithm [108] to identify the consensus sequence from the trained pHMM, which translates into the corrected assembly sequence. Apollo uses a slightly modified design of pHMMs to avoid certain limitations associated with traditional pHMMs when generating the consensus sequences [92, 93]. The modified design avoids loops in the insertion states and uses transitions to account for deletions instead of deletion states. These modifications allow the pHMM-based error correction applications [23, 24, 27] to construct more accurate consensus sequences from pHMMs.

Protein Family Search. Classifying protein sequences into families is widely used to analyze the potential functions of the proteins of interest [109, 110, 111, 112, 113, 114]. The protein family search finds the family of the protein sequence in existing protein databases. A pHMM usually represents one protein family in the database to avoid searching for many individual sequences. The protein sequence can then be assigned to a protein family based on the similarity score of the protein when compared to a pHMM in a database. This approach is used to search protein sequences in the Pfam database [115], where the HMMER [33] software suite is used to build HMMs and assign query sequences to the best fitting Pfam family. Similar to the Pfam database, HMMER’s protein family search tool is integrated into the European Bioinformatics Institute (EBI) website as a web tool. The same approach is also used in several other important applications, such as classifying many genomic sequences into potential viral families [116].

Multiple Sequence Alignment. MSA detects the differences between several biological sequences. Dynamic programming algorithms can optimally find differences between genomic sequences, but the complexity of these algorithms increases drastically with the number of sequences [117, 118]. To mitigate these computational problems, heuristic algorithms are used to obtain an approximate yet computationally efficient solution for multiple alignments of genomic sequences. pHMM-based approaches provide an efficient solution for MSA [119]. The pHMM approaches, such as hmmalign [33], assign likelihoods to all possible combinations of differences between sequences to calculate the pairwise similarity scores using forward and backward calculations or other optimization methods (e.g., particle swarm optimization [120]). pHMM-based MSA approaches are mainly useful to avoid making redundant comparisons as a sequence can be compared to a pHMM graph, similar to the protein family search.

3.1 Sources of Inefficiencies

To identify and understand the performance overheads of state-of-the-art pHMM-based applications, we thoroughly analyze existing tools for the three use cases of pHMM: (1) Apollo [24] for error correction, (2) hmmsearch [33] for protein family search, and (3) hmmalign [33] for MSA. We make six key observations based on our profiling with Intel VTune [121] and gprof [122].

Observation 1: The Baum-Welch Algorithm causes significant computational overhead. Figure 2 shows the percentage of the execution time of all three steps in the Baum-Welch algorithm for the three bioinformatics applications. We find that the Baum-Welch algorithm causes significant performance overhead for all three applications as the algorithm constitutes from 45.76% up to \(98.57\%\) of the total CPU execution time. Our profiling shows that these applications are mainly compute-bound. Forward and Backward calculations are the common steps in all three applications, whereas the Parameter Updates step is executed only for error correction. This is because protein family search and MSA use the Forward and Backward calculations mainly for scoring between a sequence and a pHMM graph as part of inference. We do not include the cost of training for these applications as it is executed only once or a few times, such that the cost of training is insignificant compared to the frequently executed inference. However, the nature of error correction requires frequently performing both training and inference for every input sequence such that the cost of training is substantial for this application. As a result, accelerating the entire Baum-Welch algorithm is key for improving the end-to-end performance of the applications.

Fig. 2.

View Figure

Observation 2: SIMD-based tools on CPU and GPUs provide suboptimal vectorization. The Baum-Welch algorithm involves frequent floating-point multiplications and additions. To resolve performance issues, several CPU-based tools (e.g., hmmalign) use SIMD instructions. However, these tools exhibit poor SIMD utilization due to inadequate port utilization and low vector capacity usage (below 50%). This suggests that CPU-based optimizations for floating-point operations, such as SIMD instructions, provide limited computational benefits for the Baum-Welch algorithm. We further investigate if the SIMD utilization in GPUs exhibits similar low utilization. To observe this, we profile our GPU work, ApHMM-GPU, to execute the two main kernels in the application: Forward and Backward calculations. We observe that the Forward calculation suffers from low SIMD utilization (i.e., percentage of active threads per warp) of around 50%, while the SIMD utilization of Backward calculation is usually close to 100%. The GPU implementation iterates over all the states that have a connection to the state that the thread is working on. However, the number of states to iterate can substantially be different per thread during Forward calculation as insertion and match states have a largely different number of incoming states, which is not the case in Backward calculation. This imbalance causes high warp divergence during Forward calculation, reducing the SIMD utilization.

Observation 3: A significant portion of floating-point operations is redundant. We observe that the same multiplications are repeatedly executed in the training step, because certain floating-point values associated with transition and emission probabilities are mainly constant during training in error correction. Our profiling analysis with VTune shows that these redundant computations constitute around \(22.7\%\) of the overall execution time when using the Baum-Welch algorithm for training in error correction.

Observation 4: Filtering the states is costly despite its advantages. The Baum-Welch algorithm requires performing many operations for a large number of states. These operations are repeated in many iterations, and the number of states can grow with each iteration. There are several approaches to keep the state space (i.e., number of states) near-constant to improve the performance or the space efficiency of the Baum-Welch algorithm [24, 27, 123, 124, 125, 126, 127]. A simple approach is to pick the best-n states that provide the highest scores in each iteration while the rest of the states are ignored in the next iteration, known as filtering [27]. Figure 3 shows the relation between the filter size (i.e., the number of states picked as best-n states), the overall runtime of the Baum-Welch algorithm, and its corresponding accuracy. Although the filtering approach is useful for reducing the runtime without significantly degrading the overall accuracy of the Baum-Welch algorithm, such an approach requires extra computations (e.g., sorting) to pick the best-n states. We find that such a filtering approach incurs substantial performance costs by constituting around \(8.5\%\) of the overall execution time in the training step.

Fig. 3.

View Figure

Observation 5: HMM accelerators are suboptimal for accelerating pHMMs. Generic HMMs do not require constraints on the connection between states (i.e., transitions) and the number of states. pHMMs are a special case for HMMs where transitions are predefined, and the number of states is determined based on the sequence that a pHMM graph represents. These design choices in HMMs and pHMMs affect the data dependency patterns when executing the Baum-Welch Algorithm. Figure 4 shows an example of the data dependency patterns in pHMMs and HMMs when executing the Baum-Welch algorithm. We observe that although HMMs and pHMMs provide similar temporal locality (e.g., only the values from the previous iteration are used), pHMMs provide better spatial locality with their constrained design. This observation suggests that HMM-based accelerators cannot fully exploit the spatial locality that pHMMs provide as they are oblivious to the design of pHMMs.

Fig. 4.

View Figure

Observation 6: Non-genomics pHMM-based applications also suffer from the computational overhead of the Baum-Welch algorithm. Among many non-genomics pHMM-based implementations [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57], we analyze the available CPU implementation of a recent pattern-matching application that uses pHMMs [52]. Our initial analysis shows that almost the entire execution time (\(98\%\)) of this application is spent on the Forward calculation, and it takes significantly longer to execute a relatively small dataset compared to the bioinformatics applications.

Many applications use either the entire or parts of the Baum-Welch algorithm for training the probabilities of HMMs and pHMMs [23, 24, 25, 27, 33, 37, 38, 40, 41, 42, 46, 47, 50, 51, 52, 53, 55, 128]. However, due to computational inefficiencies, the Baum-Welch algorithm can result in significant performance overheads on these applications. Solving the inefficiencies in the Baum-Welch algorithm is mainly important for services that frequently use these applications, such as the EBI website using HMMER for searching protein sequences in protein databases [129]. Based on the latest report in 2018, there have been more than 28 million HMMER queries on the EBI website within two years (2016–2017) [130]. These queries execute parts of the Baum-Welch algorithm more than 38,000 times daily. Such frequent usage leads to significant waste in compute cycles and energy due to the inefficiencies in the Baum-Welch algorithm.

While the Baum-Welch algorithm is computationally intensive and can consume a significant portion of the runtime and energy in various applications, these applications are often run multiple times as part of routine analyses or when new data becomes available. For error correction, the assembly of a particular genome can be reconstructed and corrected multiple times if additional sequencing data for the genome becomes available or if the tools used in the assembly construction pipeline are updated or replaced. For the protein family search and the multiple sequence alignment, protein sequencing data is frequently used multiple times due to regular updates in databases like the Pfam database [115, 131]. These updates can generate new insights [132], such as more accurate reannotation of genes in assemblies [133]. This frequent use of sequenced data can make using the Baum-Welch algorithm a time and energy-consuming process in the overall sequencing data analysis pipeline. Improving the efficiency of the Baum-Welch algorithm can significantly reduce both the compute cycles and energy consumption, especially in use cases where sequencing data is used multiple times.

3.2 Goal

Based on our observations, we find that we need to have a specialized, flexible, high-performant, and energy-efficient design to ➊ support different pHMM designs with specialized compute units for each step in the Baum-Welch algorithm, ➋ eliminate redundant operations by enabling efficient reuse of the common multiplication products, ➌ exploit spatiotemporal locality in on-chip memory, and ➍ perform efficient filtering. Such a design has the potential to significantly reduce the computational and energy overhead of the applications that use the Baum-Welch algorithm in pHMMs. Unfortunately, software- or hardware-only solutions cannot solve these inefficiencies easily. There is a pressing need to develop a hardware-software co-designed and flexible acceleration framework for several pHMM-based applications that use the Baum-Welch algorithm.

In this work, our goal is to reduce computational and energy overheads of pHMMs-based applications that use the Baum-Welch algorithm with a flexible, high-performance, energy-efficient hardware-software co-designed acceleration framework. To this end, we propose ApHMM, the first highly flexible, high-performance, and energy-efficient accelerator that can support different pHMM designs to accelerate wide-range pHMM-based applications.

4.1 Microarchitecture Overview

ApHMM provides a flexible, high-performant, and energy-efficient hardware-software co-designed acceleration framework for calculating each step in the Baum-Welch algorithm. Figure 5 shows the main flow of ApHMM when executing the Baum-Welch algorithm for pHMMs. To exploit the massive parallelism that DNA and protein sequences provide, ApHMM processes many sequences in parallel using multiple copies of hardware units called ApHMM Cores. Each ApHMM core aims to accelerate the Baum-Welch algorithm for pHMMs. An ApHMM core contains two main blocks: (1) Control Block and (2) Compute Block. Control Block provides efficient on- and off-chip synchronization and communication with CPU, DRAM, and L2/L1 cache. Compute Block efficiently and effectively performs each step in the Baum-Welch algorithm: (1) Forward calculation, (2) Backward calculation, and (3) Parameter Updates with respect to their corresponding equations in Section 2.2.

Fig. 5.

View Figure

ApHMM starts when the CPU loads necessary data to memory and sends the parameters to ApHMM ➊. ApHMM uses the parameters to decide on the pHMM design (i.e., either traditional pHMM design or modified design for error correction) and steps to execute in the Baum-Welch algorithm. The parameters related to design are sent to Compute Block so that each Compute Block can efficiently make proper state connections ➋. For each character in the input sequence that we aim to calculate the similarity score, Compute Block performs (1) Forward, (2) Backward, (3) and Parameter Updates steps. ApHMM enables disabling the calculation of Backward and Parameter Updates steps if they are not needed for an application. ApHMM iterates over the entire input sequence to fully perform the Forward calculation with respect to Equation (1)➌. ApHMM then re-iterates each character on the input sequence character-by-character to perform the Backward calculations for each timestamp t with respect to Equation (2) (i.e., step-by-step) . ApHMM updates emission and transition probabilities as the Backward values are calculated in each timestamp.

4.2 Control Block

Control Block is responsible for managing the input and output flow of the compute section efficiently and correctly by issuing both memory requests and proper commands to Compute Block to configure for the next set of operations (e.g., the forward calculation for the next character of the sequence S). Figure 5 shows three main units in Control Block: (1) Parameters, (2) Data Control, and (3) Histogram Filter.

Parameters . Control Block contains the parameters of pHMM and the Baum-Welch algorithm. These parameters define (1) pHMM design (i.e., either the traditional design or modified design for error correction) and (2) steps to execute in the Baum-Welch algorithm as ApHMM allows disabling the calculation of Backward or Parameter Updates steps.

Data Control . To ensure the correct, efficient, and synchronized data flow, ApHMM uses Data Control to (1) arbitrate among the read and write clients and (2) pipeline the read and write requests to the memory and other units in the accelerator (e.g., Histogram Filter). Data control is the main memory management unit for issuing a read request to L1 cache to obtain (1) each input sequence S, (2) corresponding pHMM graph (i.e., \(G(V,A)\)), (3) corresponding parameters and coefficients from the previous timestamp (e.g., Forward coefficients from timestamp \(t-1\) as shown in Equation (1)). Data Control collects and controls the write requests from various clients to ensure data is synchronized.

Histogram Filter . The filtering approach is beneficial for eliminating unnecessary (i.e., negligible) states from Forward and Backward calculations without significantly compromising accuracy (Section 3). The challenge in implementing a straightforward filtering mechanism lies in performing sorting in hardware, which is difficult to achieve efficiently. Our key idea is to replace the sorting mechanism with a histogram-based filter, allowing values to be placed into different bins based on their Forward or Backward values. This offers quick and approximate identification of necessary states (i.e., states with the best values until the filter is full) based on their bin locations. To enable such a binning mechanism, we employ a flexible histogram-based filtering mechanism in the ApHMM on-chip memory.

Figure 6(a) shows the overall structure of our Histogram Filter. The filter categorizes states into bins based on their Forward or Backward values at the current execution timestamp in three steps. First, Histogram Filter divides the \([0, 1]\) range into n bins, with each bin corresponding to a specific range of Forward or Backward values. The range for each bin is \(1/n\). We empirically chose 16 bins, ensuring a range of \(1/16 = 0.0625\), to maintain the same accuracy when the filter size is 500 (Figure 3). For simplicity, we use 0.06 as the range value in Figure 6(a), with the maximum value in each bin’s range displayed under Max. Value.

Fig. 6.

View Figure

Second, the Histogram Filter assigns addresses to states such that all states with Forward or Backward values within the same range fall into the same memory block. This addressing mechanism employs a base and offset strategy, where the base represents the start of the memory block for a specific range of values, and the offset is the pointer to the next free memory region within the memory block. This base and offset strategy allows ApHMM to discard unnecessary states efficiently, as their addresses are known without sorting.

Third, to identify the addresses of unnecessary states, the Histogram Filter accumulates the count of states in each bin, starting with the bin with the largest Max. Value (i.e., 1.00). When the overall state count exceeds the filter size (e.g., 500), the remaining bins are assumed to contain only unnecessary states. The Histogram Filter can find all the necessary states that a filtering technique with a sorting mechanism finds, albeit with the cost of including states beyond the predetermined filter size, as the accumulated state count in the last bin can exceed the filter count. While it is possible to perform additional computations in the last bin to prevent exceeding the filter size, we leave such optimization for future work.

To build a flexible framework for various applications, The microarchitecture is configurable to vary the number of bins (n) based on the application and the average sequence length. We recommend conducting an empirical analysis before determining this range for a particular application, as it may vary significantly.

ApHMM offers the option to disable the filtering mechanism if the application does not necessitate a filter operation for more optimal computations. Figure 6(b) shows the normalized runtime of ApHMM with and without using a filter for varying sequence lengths. We observe that the performance significantly improves when the filtering mechanism is enabled, especially for longer sequences. This can be primarily attributed to the fact that the number of states requiring processing at the subsequent timestamp can exponentially increase, as each state typically has more than one transition, potentially leading to an exponential increase in states at each subsequent timestamp. As the sequence length grows, such an exponential increase can adversely affect the application, which can be significantly mitigated without compromising accuracy through a filtering approach.

4.3 Compute Block

Figure 7 shows the overall structure of a Compute Block, which is responsible for performing core compute operations of each step in the Baum-Welch algorithm (Figure 5) based on the configuration set by the Control Block via Index Control➊. A Compute Block contains two major units: (1) a unit for calculating Forward (Equation (1)) and Backward (Equation (2)) values➋ and updating transition probabilities (Equation (3)) , and (2) a unit for updating the emission probabilities (Equation (4)) . Each unit performs the corresponding calculations in the Baum-Welch algorithm.

Fig. 7.

View Figure

Forward and Backward Calculations. Our goal is to calculate the forward and backward values for all states in a pHMM graph \(G(V,A)\), as shown in Equations (1) and (2), respectively. To calculate the Forward or Backward value of a state i at a timestamp t, ApHMM uses Processing Engines (PEs). Since pHMMs may require processing hundreds to thousands of states to process at a time, ApHMM includes many PEs and groups them into PE Groups. Each PE is responsible for calculating the Forward and Backward values of a state \(v_i\) per timestamp t. Our key challenge is to balance the utilization of the compute units with available memory bandwidth. We discuss this trade-off between the number of PEs and memory bandwidth in Section 4.4. To efficiently calculate the Forward and Backward values, PE performs two main operations.

First, PE uses the parallel four lanes in Dot Product Tree and Accumulator to perform multiple multiply and accumulation operations in parallel, where the final summation is calculated in the Reduction Tree. This design enables efficient multiplication and summation of values from previous timestamps (i.e., \(F_{t-1}(j)\) or \(B_{t+1}(j)\)). Second, to avoid redundant multiplications of transition and emission probabilities, the key idea in PEs is to efficiently enable the reuse of the products of these common multiplications. To achieve this, our key mechanism stores these common products in lookup tables (LUTs) in each PE while enabling efficient retrievals of the common products. We store these products as these values can be preset (i.e., fixed) before the training step starts and frequently used during training while causing high computational overheads.

Our key challenge is to design minimal but effective LUTs to avoid area and energy overheads associated with LUTs without compromising the computational efficiency LUTs provide. To this end, we analyze error correction, protein family search, and multiple sequence alignment implementations. We observe that (1) redundant multiplications are frequent only during training and (2) the alphabet size of the biological sequence significantly determines the number of common products (i.e., 4 in DNA and 20 in proteins). Since error correction is mainly bottlenecked during the training step, we focus on the DNA alphabet and the pHMM design that error correction uses. We identify that each state uses (1) at most 4 different emission probabilities (i.e., DNA letters) and (2) on average 7 different transitions. This results in 28 different combinations of emission and transition probabilities. To enable slightly better flexibility, we assume 9 different transitions and include 36 entries in LUTs.

The key benefit is LUTs provide ApHMM with a bandwidth reduction of up to \(66\%\) per PE while avoiding redundant computations. ApHMM is flexible such that it enables disabling the use of LUTs and instead performing the actual multiplication of transition and emission probabilities (TE MUL unit in Figure 7).

Updating the Transition Probabilities. Our goal is to update the transition probabilities of all the states, as shown in Equation (3). To achieve this, we design the Update Transition(UT) compute unit and tightly couple it with PEs, as shown in Figure 7. Each UT efficiently calculates the denominator and numerator in Equation (3) for a state \(v_i\). UTs include three key mechanisms.

First, to enable efficient broadcasting of common values between the Backward calculation and Parameter Updates steps, ApHMM connects PEs with UTs for updating transitions. Each PE in a PE Group is broadcasted with the same previously calculated \(F_{t}(i)\) or \(B_{t+1}(j)\) values from the previous timestamp for calculating the \(F_{t+1}(j)\) or \(B_{t}(i)\) values, respectively. The incoming red arrows in Figure 7 show the flow of these Forward and Backward values in PEs and UTs. The calculation of \(F_{t+1}(j)\) involves a summation over all states i as shown in Equation (1). The \(F_{t}(i)\) term is common to the calculation of \(F_{t+1}(j)\) for all states j and hence can be broadcasted. Similarly, the calculation of \(B_{t}(i)\) involves a summation over all states j (Equation (2)). The \(B_{t+1}(j)\) term is common to the calculation of \(B_{t}(i)\) for all states i and hence can be broadcasted. This key design choice exploits the broadcast opportunities available within the common multiplications in the Baum-Welch equations.

Second, ApHMM cores are designed to directly consume the broadcasted Backward values when updating the emission and transition probabilities to reduce the bandwidth and storage requirements. We exploit the broadcasting opportunities, because we observe that Backward values do not need to be fully computed, and they can be directly consumed when updating the transitions and emission probabilities while the Backward values are broadcasted in the current timestamp. ApHMM updates emission and transition probabilities step-by-step as Backward values are calculated, a hardware-software optimization we call the partial compute approach. It is worth noting that ApHMM fully computes and stores the Forward values before updating the emission and transition probabilities. The key benefits of combining broadcasting with the partial compute approach are (1) decoupling hardware scaling from bandwidth requirements and (2) reducing the bandwidth requirement by \(4\times\) (i.e., 32 bits/cycle instead of 128 bits/cycle).

Third, to exploit the spatiotemporal locality in pHMMs, we utilize on-chip memory in UTs with memoization techniques that allow us to store the recent transition calculations. We observe from Equation (3) that transition update is calculated using the values of states connected to each other. Since the connections are predefined and provide spatial locality (Figure 4), our key idea is to memoize the calculation of all the numerators from the same i to different states by storing these numerators in the same memory space. This enables us to process the same state i in different timestamps within the same PE to reduce the data movement overhead within ApHMM. To this end, we use an 8 KB on-chip memory (Transition Scratchpad) to store and reuse the result of the numerator of Equation (3). Since we store the numerators that contribute to all the transitions of a state i within the same memory space, we perform the final division in Equation (3) by using the values in the Transition Scratchpad. We use an 8 KB memory as this enables us to store 256 different numerators from any state i to any other state j. We observe that pHMMs have 3–12 distinct transitions per state. Thus, 8 KB storage enables us to operate on at least 20 different states within the same PE. The memoization technique allows (1) skipping redundant data movement and (2) reducing the bandwidth requirement by \(2 \times\) per UT.

Updating the Emission Probabilities. Our goal is to update the emission probabilities of all the states, as shown in Equation (4). To achieve this, we use the Update Emission (UE) unit, as shown in Figure 7, which includes three smaller units: (1) Calculate Emission Numerator, (2) Calculate Emission Denominator, and (3) Division & Update Emission. UE performs the numerator and denominator computations in parallel as they are independent of each other, which includes a summation of the products \(F_{t}(i)B_t(i)\). These \(F_{t}(i)\) and \(B_t(i)\) values are used to update both the transition and emission probabilities, as shown in Equation (3). To reduce redundant computations, our key design choice is to use the \(F_{t}(i)\) and \(B_t(i)\) values as broadcasted in the transition update step since these values are also used for updating the emission probabilities. Thus, we broadcast these values to UEs through Write Selectors, as shown in Figure 7.

An ApHMM core writes and reads both the numerator and denominator values to L1 cache to update the emission probabilities. The results of the division operations and the posterior emission probabilities (i.e., \(e^{*}_{X}(v_{i})\) in Equation (4)) are written back to L1 cache after processing each read sequence S. If we assume that the number of characters in an alphabet \(\Sigma\) is \(n_{\Sigma }\) (e.g., \(n_{\Sigma } = 4\) for DNA letters), then ApHMM stores \(n_{\Sigma }\) many different numerators for each state of the graph as emission probability may differ per character for each state. Our microarchitecture design is flexible such that it allows defining \(n_{\Sigma }\) as a parameter.

4.4 Hardware Configuration Choice

Our goal is to identify the ideal number of memory ports and PEs for better scaling ApHMM with many cores. We identify the number of memory ports and their dependency on the hardware scaling in four steps. First, ApHMM requires one input memory port for reading the input sequence to update the probabilities in a pHMM graph. Second, updating the transition probabilities requires 3 memory ports: (1) reading the forward value from L1, (2) reading the transition, and (3) emission probabilities if using the LUTs is disabled (Section 4.3). Since these ports are shared across each PE, the number of PEs and memory bandwidth per port determines the utilization of these memory ports. Third, ApHMM requires 4 memory ports to update the emission probabilities for (1) calculating the numerator and (2) denominator in Equation (4), (3) reading the forward from Write Selectors, and (4) writing the output. These memory ports are independent of the impact of the number of PEs in a single ApHMM core. Fourth, ApHMM does not require additional memory ports for each step in the Baum-Welch algorithm due to the broadcasting feature of ApHMM (Section 4.3). Instead, computing these steps depends on (1) memory bandwidth per port, which determines the number of multiplications and accumulations in parallel in a PE, and (2) the number of processing engines (PEs). We conclude that the overall requirement for a ApHMM core is eight memory ports with the same bandwidth per port.

Figure 8 shows the scaling capabilities of ApHMM with the number of PEs and sequence length to decide (1) the overall number of PEs and (2) the longest chunk size for the best acceleration. First, to decide the overall number of PEs to use in ApHMM, in Figure 8(a), we show the acceleration speedup while scaling ApHMM with the number of PEs and bandwidth per memory port, where we keep the number of memory ports fixed to 8. We observe that a linear trend of increase in acceleration is possible until the number of PEs reaches 64, where the rate of acceleration starts declining. We explore the reason for such a trend in Figure 8(b). We find that the acceleration on the transition step starts settling down as the number of PEs grows due to memory port limitation that reduces parallel data read from memory per PE, eventually resulting in the underutilization of resources. Second, We conclude that the acceleration trend we observe in Figure 8(a) is mainly due to the scaling impact on the forward and backward calculation when the number of PEs is greater than 64 where 8 memory ports start becoming the bottleneck.

Fig. 8.

View Figure

In our design, the choice of memory bandwidth influences the number of PE Groups and PEs, given a constant number of memory ports. While our hardware can scale to accommodate higher bandwidths, we opt for a 16-bytes/cycle bandwidth. This design choice aligns with the 128-bit line size of our L1 cache, allowing us to operate on four single-precision floating-point values (32-bit) across 4 PEs simultaneously. To fully utilize all 64 PEs, as discussed earlier, we employ 16 PE Groups (64 PEs = 4 PEs × 16 PE Groups).

Second, to identify the optimal chunk size (i.e., sequence length) that ensures a near-linear increase in execution time with increasing sequence length, we examine the execution time of the Baum-Welch algorithm for chunk sizes of 150, 650, and 1,000 bases, as shown in Figure 8(c). We observe a linear increase in execution time with chunk sizes up to approximately 650 bases. Beyond this point, the execution time begins to increase non-linearly. This non-linear scaling is primarily due to the increased cache space requirements for storing certain parameters (e.g., emission values), as shown in Supplemental Figure S1. This increased cache pressure leads to more data accesses from the upper levels of the memory hierarchy. ApHMM can maintain a linear trend in execution time for longer sequences by either increasing the L1 and L2 cache capacities or utilizing higher-bandwidth memories to mitigate the data movement overheads. We provide further details regarding the data distribution and memory layout in Supplemental Section S2.

We conclude that the memory ports and chunk size primarily constrain the acceleration speedup of ApHMM, as the PEs start to be underutilized due to increased data movement overheads. To further enhance the acceleration with ApHMM, optimizing the utilization of PEs by minimizing these overheads is crucial.

Number of ApHMM Cores. We show our methodology for choosing the ideal number of cores in ApHMM for accelerating the applications. Figure 9 shows the speedup of three bioinformatics applications when using single, 2, 4, and 8 cores in ApHMM. We divide the entire execution time of the applications into three stages: (1) the CPU execution of the application that does not use the Baum-Welch execution, (2) the Baum-Welch execution accelerated using ApHMM, and (3) and the overhead mainly caused due to data movements. Our analysis incorporates the estimated off- and on-chip data movement overhead. We observe that using four cores in ApHMM provides the best speedup overall. This is because the applications provide smaller rooms for acceleration for two reasons. First, the remaining CPU part of the application becomes the bottleneck in the entire execution of the application due to the significant acceleration of the Baum-Welch execution using ApHMM. Second, the data movement overhead starts causing more latency than the benefits of further accelerating the Baum-Welch algorithm by increasing the number of cores. This suggests ApHMM is bounded by the data movement overhead when scaling it to a larger number of cores, and there is still room for improving the performance of ApHMM by placing ApHMM inside or near the memory (e.g., high-bandwidth memories) to reduce the data movement overheads that limit scaling ApHMM to many cores. Based on our observations, we use a four-core ApHMM to achieve the best overall performance (see Supplemental Section S3 for the execution flow of the system with multiple cores in ApHMM).

Fig. 9.

View Figure

5.1 Evaluation Methodology

We use the configurations shown in Table 1 to implement the ApHMM design described in Section 4 in SystemVerilog. We carry out synthesis using Synopsys Design Compiler [134] in a typical 28 nm process technology node at 1 GHz clock frequency with tightly integrated on-chip memory (1 GHz) to extract the logic area and power numbers. We develop an analytical model to extract performance and area numbers for a scale-up configuration of ApHMM. We use four ApHMM cores in our evaluation (Section 4.4). We account for an additional 5% of cycles to compensate for the arbitration across memory ports. These extra cycles estimate the cycles for synchronously loading data from DRAM to L2 memory of a single ApHMM core and asynchronous accesses to DRAM when more data needs to be from DRAM for a core (e.g., Forward calculation may not fit the L2 memory).

We use the CUDA library [135] (version 11.6) to provide a GPU implementation of the software optimizations described in Section 4 for executing the Baum-Welch algorithm. Our GPU implementation, ApHMM-GPU, uses the pHMM design designed for error correction, implements LUTs (Section 4.3) as shared memory, and uses buffers to arbitrate between current and previous Forward/Backward calculations to reflect the software optimizations of ApHMM in GPUs. We integrate our GPU implementation with a pHMM-based error correction tool, Apollo [24], to evaluate the GPU implementation. Our GPU implementation is the first GPU implementation of the Baum-Welch algorithm for profile Hidden Markov models.

We use gprof [122] to profile the baseline CPU implementations of the use cases on the AMD EPYC 7742 processor (2.26 GHz, 7 nm process) with single- and multi-threaded settings. We use the CUDA library and nvidia-smi to capture the runtime and power usage of ApHMM-GPU on NVIDIA A100 and NVIDIA Titan V GPUs, respectively.

We compare ApHMM with the CPU, GPU, and FPGA implementations of the Baum-Welch algorithm and use cases in terms of execution time and energy consumption. To evaluate the Baum-Welch algorithm, we execute the algorithm in Apollo [24] and calculate the average execution time and energy consumption of a single execution of the Baum-Welch algorithm. To evaluate the end-to-end execution time and energy consumption of error correction, protein family search, and multiple sequence alignment, we use Apollo [24], hmmsearch [33], and hmmalign [33]. We replace their implementation of the Baum-Welch algorithm with ApHMM when collecting the results of the end-to-end executions of the use cases accelerated using ApHMM. When available, we compare the use cases that we accelerate using ApHMM to their corresponding CPU, GPU, and FPGA implementations. For the GPU implementations, we use both ApHMM-GPU and HMM_cuda [96]. For the FPGA implementation, we use the FPGA Divide and Conquer (D&C) accelerator proposed for the Baum-Welch algorithm [95]. When evaluating the FPGA accelerator, we ignore the data movement overhead and estimate the acceleration based on the speedup provided by the earlier work. We acknowledge that the performance and energy comparisons can be attributed to both platform differences and architectural optimizations, especially when comparing ApHMM with the FPGA accelerator. Although our evaluations lack comparisons in the equivalent platforms for FPGAs, we still believe that our evaluations provide valuable insights regarding the benefits of our ASIC implementation compared to the FPGA work.

In terms of accuracy, we ensure the accuracy of our results by faithfully implementing all the equations of the Baum-Welch algorithm and rigorously testing their output during our ASIC design. The only exception is the Histogram Filter, where we introduce a binning approach to include all the states a sorting-based software implementation would include, ensuring at least the same minimum accuracy as the original software implementation. Our accuracy evaluation shows that the histogram filter approach usually leads to better accuracy than the sorting approach, with a minimal accuracy difference within a \(+/-\)0.2% range. To reproduce the output for comparison purposes, we provide the source code of our software optimizations in the GPU implementation of ApHMM (ApHMM-GPU).

Data Set. To evaluate the error correction use case, we prepare the input data that Apollo requires: (1) assembly and (2) read mapping to the assembly. To construct the assembly and map reads to the assembly, we use reads from a real sample that includes overall 163,482 reads of Escherichia coli(E. coli) genome sequenced using PacBio sequencing technology with the average read length of 5,128 bases. The accession code of this sample is SAMN06173305. Of 163,482 reads, we randomly select 10,000 sequencing reads for our evaluation. We use minimap2 [136] and miniasm [137] to (1) find overlapping reads and (2) construct the assembly from these overlapping reads, respectively. To find the read mappings to the assembly, we use minimap2 to map the same reads to the assembly that we generate using these reads. We provide these inputs to Apollo for correcting errors in the assembly we construct.

To evaluate the protein family search, we use the protein sequences from a commonly studied protein family, Mitochondrial carrier (PF00153), which includes 214,393 sequences with an average length of 94.2. We use these sequences to search for similar protein families from the entire Pfam database [115] that includes 19,632 pHMMs. To achieve this, the hmmsearch [33] tool performs the Forward and Backward calculations to find similarities between pHMMs and sequences.

To evaluate the multiple sequence alignment, we use 1,140,478 protein sequences from protein families Mitochondrial carrier (PF00153), Zinc finger (PF00096), bacterial binding protein-dependent transport systems (PF00528), and ATP-binding cassette transporter (PF00005). We align these sequences to the pHMM graph of the Mitochondrial carrier protein family. To achieve this, the hmmalign [33] tool performs the Forward and Backward calculations to find similarities between a single pHMM graph and sequences.

Data Set Justification. In our study, we carefully chose our datasets for overhead analysis and evaluation. We believe these datasets are comprehensive and relevant to guide our ASIC design and to evaluate ApHMM with other systems for several reasons. First, our datasets cover various use cases with various sequence lengths (i.e., an average read length of 5,168 and an average protein sequence length of 94.2) and alphabet sizes (4 in DNA and 20 in proteins). This diversity ensures that our results are not skewed toward a specific use case or dataset.

Second, for error correction, we use a real-world sample of the E. coli genome, a commonly studied bacterial genome. The overall length of randomly selected 10,000 E. coli reads is around 50,000,000 bases (the average length of a single read is 5,168). This ensures that these reads cover the entire E. coli genome around 10 times (i.e., \(10\times\) depth of coverage), which ensures that the Baum-Welch algorithm is executed by performing error correction on the entire genome multiple times without focusing on the specific regions of the genome to avoid potential bias that can be caused on particular regions. For the multiple sequence alignment and the protein family search, we use the most commonly studied protein families as these protein families are among the top 20 families with the largest number of protein sequence alignments,¹ ensuring the relevance and applicability of our work. The bottleneck analysis was conducted on a subset of these datasets, demonstrating that our design is not overfitting to a specific dataset.

Third, the Baum-Welch algorithm operates on a sub-region of the pHMM graph, the size of which is determined by the sequence length or chunk size, whichever is shorter. Thus, the complexity of a single Baum-Welch execution on this sub-region is determined mainly based on the sequence length and the alphabet size, regardless of the overall genome size or the sequence lengths larger than the chunk size. In our case, we cover all these cases: (1) the pHMM subgraph is determined based on the sequence length (around 90 bases) as it is shorter than the chunk size (up to 1,000 bases) in the protein family search and the multiple sequence alignment (2) the length of the pHMM subgraph is determined by the chunk size in error correction as the sequence length is usually larger (around 10,000 bases) than the chunk size, and (3) different alphabet sizes in DNA and protein.

Fourth, for overhead analysis, we discuss in Section 3.1, we ensure our design is not overfitting to a specific dataset by using a subset of these datasets for each use case. The overhead was measured by taking the geometric mean across different runs to further ensure the robustness of our design. Since our ASIC design is mainly influenced based on the observations we make in our overhead analysis and to maximize the performance improvement for the applications mainly bottlenecked by the Baum-Welch algorithm (i.e., error correction), we believe the comprehensiveness of our data set choice and the overhead analysis enable us improving the robustness of our accelerator across a wide range of potentially many other use cases other than the use cases we evaluate in this work.

5.2 Area and Power

Table 2 shows the area breakup of the major modules in ApHMM. For the area overhead, we find that the UT units take up most of the total area (\(77.98\%\)). This is mainly because UTs consist of several complex units, such as a multiplexer, division pipeline, and local memory. For the power consumption, Control Block and PEs contribute to almost the entire power consumption (\(86\%\)) due to the frequent memory accesses these blocks make. Overall, aApHMM core incurs an area overhead of 6.5 mm² in 28 nm with a power cost of 0.509 W.

5.3 Accelerating the Baum-Welch Algorithm

Figure 10 shows the performance and energy improvements of ApHMM for executing the Baum-Welch algorithm. Based on these results, we make six key observations. First, we observe that ApHMM is \(15.55\times\)–\(260.03\times\), \(1.83\times\)–\(5.34\times\), and \(27.97\times\) faster than the CPU, GPU, and FPGA implementations of the Baum-Welch algorithm, respectively. Although our evaluations do not directly compare the state-of-the-art FPGA work with the potential FPGA implementation of ApHMM, we believe the performance benefits that ApHMM provides arise not only from the differences in the platform and architecture but also from the optimizations we provide, which are absent in the existing FPGA work. We believe the benefits of these optimizations on the same platform can partly be observed when comparing ApHMM-GPU with the state-of-the-art GPU accelerator. Second, ApHMM reduces the energy consumption for calculating the Baum-Welch algorithm by \(2474.09\times\) and \(896.70\times\)–\(2622.94\times\) compared to the single-threaded CPU and GPU implementations, respectively. These speedups and reduction in energy consumption show the combined benefits of our software-hardware optimizations. Third, the parameter update step is the most time-consuming step for the CPU and the GPU implementations, while ApHMM takes the most time in the forward calculation step. The reason for such a trend shift is that ApHMM reads and writes to L2 Cache and DRAM more frequently during the forward calculation than the other steps, as ApHMM requires the forward calculation step to be fully completed and stored in the memory before moving to the next steps as we explain in Section 4.3. Fourth, we observe that ApHMM-GPU performs better than HMM_cuda by \(2.02\times\) on average. HMM_cuda executes the Baum-Welch algorithm on any type of hidden Markov model without a special focus on pHMMs. As we develop our optimizations based on pHMMs, ApHMM-GPU can take advantage of these optimizations for more efficient execution. Fifth, both ApHMM-GPU and HMM_cuda provide better performance for the Forward calculation than ApHMM. We believe the GPU implementations are a better candidate for applications that execute only the Forward calculations as ApHMM targets, providing the best performance for the complete Baum-Welch algorithm. Sixth, the GPU implementations provide a limited speedup over the multi-threaded CPU implementations mainly because of frequent access to the host for synchronization and sorting (e.g., the filtering mechanism). These required accesses from GPU to host can be minimized with a specialized hardware design, as we propose in ApHMM for performing the filtering mechanism. We conclude that ApHMM provides substantial improvements, especially when we combine speedups and energy reductions for executing the complete Baum-Welch algorithm compared to the CPU and GPU implementations, which makes it a better candidate to accelerate the applications that use the Baum-Welch algorithm than the CPU, GPU, and FPGA implementations.

Fig. 10.

View Figure

Breakdown of the optimizations benefits. Table 3 shows the performance improvements that each ApHMM optimization contributes for executing the Baum-Welch algorithm given the single-core hardware configuration we discuss in Section 4.4 compared to the CPU baseline of the Baum-Welch algorithm. We estimate the speedup of Histogram Filter by eliminating the sorting mechanism from filtering while considering the overhead of redundant states included in Histogram Filter. For other optimizations, we conservatively estimate the performance speedups by considering the memory bandwidth reductions that each optimization provides, as discussed in Section 4, and the relation between acceleration speedup and the memory bandwidth requirements (Figure 8). We make five key observations. First, we find almost all optimizations contribute significantly to reducing the overall execution time of the Baum-Welch algorithm. Although Histogram Filter provides a limited speedup, this is because it constitutes around \(8.5\%\) of the overall execution time (Observation 4 in Section 3.1). Second, the tight coupling of the broadcasting and the partial compute approach provides the most significant speedups by avoiding a large number of useless data movements with significant memory bandwidth reductions. Third, the speedup from LUTs is mainly achieved by eliminating many single-precision floating-point operations, causing around \(22.7\%\) of the total execution time (Observation 3 in Section 3.1). Fourth, the speedups with the memoization technique are purely achieved by significantly reducing the data movement latency when frequently calculating the transition probabilities. Fifth, we find that the memoization and the partial compute optimizations are utilized only in the training step, and the LUTs can be useful when the alphabet size is small (e.g., 4 in DNAs) due to storage limitations, which is usually the case when the Baum-Welch algorithm is used mainly for inference with the protein sequencing data. Although these benefits cannot be fully utilized in such cases, the remaining optimizations still provide a significant speedup up to \(3.63\times\). We conclude that our optimizations achieve significant speedups for various use cases, from training with DNA sequencing data to inferring with protein sequencing data, allowing the acceleration of many applications that use the Baum-Welch algorithm with pHMMs.

5.4 Use Case 1: Error Correction

Figures 11 and 10 show the end-to-end execution time and energy reduction results for error correction, respectively. We make four key observations. First, we observe that ApHMM is \(2.66\times\)–\(59.94\times\), \(1.29\times\)–\(2.09\times\), and \(7.21\times\) faster than the CPU, GPU, and FPGA implementations of Apollo, respectively. Second, ApHMM reduces the energy consumption by \(64.24\times\) and \(71.28\times\)–\(115.46\times\) compared to the single-threaded CPU and GPU implementations. These two observations are in line with the observations we make in Section 5.3 as well as the motivation results we describe in Section 3: Apollo is mainly bounded by the Baum-Welch algorithm, and ApHMM accelerates the Baum-Welch algorithm significantly, providing significant performance improvements and energy reductions for error correction. We conclude that ApHMM significantly improves the energy efficiency and performance of the error correction mainly because the Baum-Welch algorithm constitutes a large portion of the entire use case.

Fig. 11.

View Figure

5.5 Use Case 2: Protein Family Search

Our goal is to evaluate the performance and energy consumption of ApHMM for the protein family search use case, as shown in Figures 11 and 10, respectively. We make three key observations. First, we observe that ApHMM provides speedup by \(1.61\times\)–\(1.75\times\), and \(1.03\times\) compared to the CPU and FPGA implementations. Second, ApHMM is \(1.75\times\) more energy efficient than the single-threaded CPU implementation. The speedup ratio that ApHMM provides is lower in the protein family search than error correction, because (1) ApHMM accelerates a smaller portion of the protein family search (\(45.76\%\)) than error correction (\(98.57\%\)), and (2) the protein alphabet size (20) is much larger than the DNA alphabet size (4), which increases the DRAM access overhead of ApHMM by \(12.5\%\). Due to the smaller portion that ApHMM accelerates and increased memory accesses, it is expected that ApHMM provides lower performance improvements and energy reductions compared to the error correction use case. Third, ApHMM can provide better speedup compared to the multi-threaded CPU as a large portion of the parts that ApHMM does not accelerate can still be executed in parallel using the same amount of threads, as shown in Figure 11. We conclude that ApHMM improves the performance and energy efficiency for the protein family search, while there is a smaller room for acceleration compared to the error correction.

5.6 Use Case 3: Multiple Sequence Alignment

Our goal is to evaluate the ApHMM’s end-to-end performance and energy consumption for the MSA, as shown in Figures 11 and 10, respectively. We make three key observations. First, we observe that ApHMM performs \(1.95\times\) and \(1.03\times\) better than the CPU and FPGA implementations, while ApHMM is \(1.96\times\) more energy efficient than the CPU implementation of MSA. We note that the hmmalign tool does not provide the multi-threaded CPU implementation for MSA. ApHMM provides better speedup for MSA than the protein family search, because MSA performs more forward and backward calculations (\(51.44\%\)) than the protein search use case (\(45.76\%\)), as shown in Figure 2. Third, ApHMM provides slightly better performance than the existing FPGA accelerator (FPGA D&C) in all applications, even though we ignore the data movement overhead of FPGA D&C, which suggests that ApHMM may perform much better than FPGA D&C in real systems. We conclude that ApHMM improves the performance and energy efficiency of the MSA use case better than the protein family search.