Heterogeneous acceleration system for bwa-mem sequence alignment software and implementation method thereof

Abstract

The application provides a heterogeneous acceleration system for BWA-MEM sequence alignment software and an implementation method thereof, comprising a CPU program and an FPGA accelerator; the CPU program adopts a heterogeneous batch processing strategy, parses an input alignment file, sequentially starts the FPGA accelerator, processes FPGA accelerator output, executes a seed filtering algorithm and a seed chain generation algorithm in BWA-MEM; the heterogeneous batch processing strategy performs batch processing reconstruction on the SMEM search algorithm in the seed generation stage and the Smith-Waterman algorithm in the seed extension stage. In the case that the output result is completely same as the original BWA-MEM software, the BWA-MEM software alignment speed is greatly improved, and the sequence alignment throughput is improved.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, specifically to a heterogeneous acceleration system for BWA-MEM sequence alignment software and its implementation method. Background Technology

[0002] Gene sequencing technology can extract the complete sequence of an individual's genome for biological analysis and interpretation, identify disease-causing genes, discover potential disease risks, and enable early prevention and treatment. Next-generation sequencing technology, also known as high-throughput sequencing, can simultaneously obtain the sequences of hundreds of thousands to millions of nucleic acid molecules in a single run. Genetic data requires further analysis using alignment algorithms before it can be used. Currently, mainstream alignment algorithms include dynamic programming-based algorithms and algorithms based on auxiliary data structures. Dynamic programming-based alignment algorithms mainly include the Smith-Waterman (SW) algorithm and the Needleman-Wunsch (NW) algorithm. These algorithms calculate a score matrix between the sequence to be aligned and the reference genome to obtain their similarity. In practical applications, achieving a complete score matrix calculation requires significant computational, storage, and processing time. To address these issues, algorithms based on auxiliary data structures, such as FM-index and hash tables, compress genomic data using auxiliary data structures and employ corresponding algorithms to achieve accurate alignment.

[0003] Patent document CN114064551A (application number: CN202210046617.1) discloses a method for accelerating high-concurrency sequence alignment computation based on CPU+GPU heterogeneity, comprising the following steps: BWA-MEM algorithm code reconstruction; concurrent task processing on the CPU: completing the partitioning of the sequence set, forming multiple concurrent tasks for the first time; running the BWA-MEM algorithm after code reconstruction to complete concurrent data processing on the GPU; concurrent task processing on the GPU: for the seed set and chain generated during the sequence data comparison process, seed sets with the same or adjacent length, position, and number are partitioned into the same data block, and the chain is processed in the same way, thereby completing the partitioning of seed set and chain, forming multiple concurrent tasks for the second time.

[0004] Single alignment algorithms all have certain limitations during execution. Most alignment tools widely adopt the Seed-and-Extend paradigm, combining multiple algorithms to complete sequence alignment. The Seed-and-Extend paradigm includes two stages: seed generation and seed extension. The seed generation stage typically uses auxiliary data structures to implement an accurate alignment algorithm, finding a fragment in the genome that is identical to the sequence to be aligned on the reference genome. The seed extension stage usually uses the less precise Smith-Waterman algorithm to extend the seed, calculating the similarity between the sequence to be aligned and a fragment in the reference genome. Based on the similarity score, the sequence to be aligned is mapped to the reference genome, completing the sequence alignment and determining for information such as base variations. Among current alignment tools, BWA-MEM is the most widely used. BWA-MEM uses the SMEM search algorithm based on the FM-index data structure in the seed generation stage and the Smith-Waterman algorithm in the seed extension stage. These two algorithms are the performance bottlenecks in the complete software package's computation process.

[0005] However, the massive amounts of gene data required for comparison place higher demands on the processing speed of tools, making the acceleration of sequence alignment software a primary goal. Increasing research focuses on accelerating this process using heterogeneous computing hardware, including GPUs and FPGAs. Graphics Processing Units (GPUs), with their abundant computing core resources, are used for parallel processing of multiple sequence alignments, and this approach is widely adopted in current high-performance computing clusters for accelerating gene sequencing. However, the complex branches and dependencies inherent in sequence alignment make parallelization difficult, leading to severe thread splitting during GPU acceleration and impacting performance. FPGAs, with their pipelined architecture, can effectively handle acceleration of deeply looping serial algorithms, batch processing multiple sequences to improve throughput, making them the primary choice for current alignment algorithm acceleration hardware. However, current work using FPGAs to accelerate sequence alignment software mostly focuses on a single algorithm. A complete BWA-MEM consists of multiple algorithms, and acceleration of a single algorithm cannot work in conjunction with the software package. Furthermore, without optimized software batch processing, the data transfer latency between the FPGA and CPU is relatively high. Summary of the Invention

[0006] To address the shortcomings of existing technologies, the purpose of this invention is to provide a heterogeneous acceleration system and its implementation method for BWA-MEM sequence alignment software.

[0007] The heterogeneous acceleration system for BWA-MEM sequence alignment software provided by the present invention includes a CPU program and an FPGA accelerator.

[0008] The CPU program adopts a heterogeneous batch processing strategy to parse the input comparison file, start the FPGA accelerator in sequence, process the FPGA accelerator output, and execute the seed filtering algorithm and seed chain generation algorithm in BWA-MEM.

[0009] The heterogeneous batch processing strategy reconstructs the SMEM search algorithm in the seed generation stage and the Smith-Waterman algorithm in the seed extension stage.

[0010] The SMEM search algorithm executes this part of the program in advance. All input sequences to be compared are transmitted to the FPGA, and the accelerator calculates and obtains the super maximum exact matching seed for all sequences.

[0011] The Smith-Waterman algorithm pre-executes the bidirectional extension Smith-Waterman algorithm, batch-processes all seed chains in a sequence and the bidirectional extension of all seeds within them, and then writes the score results back to the CPU for further filtering.

[0012] Preferably, the FPGA accelerator includes an SMEM search algorithm accelerator;

[0013] The SMEM search algorithm accelerator includes a register configuration and state machine module, a forward search and backward search execution module, and a sequence and intermediate result storage and control module.

[0014] The register configuration and state machine module provides internal initialization data configuration for the accelerator, records the working state of the accelerator, and controls the accelerator's operation process through the state machine.

[0015] The forward search and backward search execution modules decouple the forward and backward search data dependencies within the software, enabling them to run in parallel. Within a single search module, the order of computation and memory access is exchanged, realizing a fully pipelined internal design structure.

[0016] The sequence and intermediate result storage and control module uses BRAM to store the bases of the sequence to be compared and the intermediate results during the search process in a dual-interval manner, and controls the writing and reading of data according to the read and write enable.

[0017] Preferably, the FPGA accelerator also includes a Smith-Waterman algorithm accelerator;

[0018] The Smith-Waterman algorithm accelerator includes a register configuration and state machine module and a systolic array module;

[0019] The pulse array module includes a PE unit array module, a pulse array controller, a data router, and a score calculator module;

[0020] The PE cell array module calculates the score of a matrix grid point within one period according to the Smith-Waterman algorithm score calculation rules.

[0021] The pulse array controller is responsible for providing enable signals to all other modules periodically and precisely.

[0022] The data router is responsible for providing different initialization data to different PE unit array modules, as well as the base data required by the PE unit array modules during operation.

[0023] The score calculator is responsible for recording the output of the score results of the PE unit array module, updating the current maximum score in each cycle, and recording the row and column where the maximum score is located.

[0024] Preferably, two hotspot functions are executed in advance: in the mem_collect_intv_fpga() function, the SMEM search algorithm accelerator on the FPGA is started to batch process all input sequences to be compared, and in the ksw_extend2_fpga() function, the Smith-Waterman algorithm accelerator on the FPGA is started to batch process the bidirectional extension of all seeds in a sequence.

[0025] Preferably, the Intel acceleration stack is used to expose the interface controller circuit of the FPGA acceleration card to the user-defined acceleration function area through bus encapsulation, the PCIe controller interface is encapsulated as the CCI-P bus protocol, and the DDR controller interface is encapsulated as the AVMM bus protocol.

[0026] An acceleration system is built based on the corresponding bus protocol and the IP modules provided by the acceleration stack. A DMA module is integrated for data transmission between the CPU and FPGA. A CCI-P bus one-to-four-converter bridge is integrated for mounting different accelerator modules. The CPU configures the internal registers of the accelerator through the CCI-P bus protocol converted from PCIe.

[0027] The acceleration system is equipped with an SMEM search algorithm accelerator and a SW algorithm accelerator. The accelerators and peripheral circuit components are clock isolated by a CCI-P asynchronous bridge and operate at an independent clock frequency. The peripheral circuit components have the same clock as the CCI-P bus protocol.

[0028] Preferably, in the SMEM search algorithm accelerator, only a specific number of sequence alignments are completed at a time. The accelerator obtains in advance the number of times the execution module needs to be called in a loop to execute all current sequences, and completes one SMEM search algorithm execution module call from the INIT state to the UPDATE state.

[0029] In the LOAD state, all the sequences to be compared that need to be processed in the current batch are moved from DDR to the storage module with the sequence to be compared. In the SMEM state, the execution module is started and the forward and backward search modules are started in sequence to calculate all the sequences to be compared in the current batch in a pipeline manner.

[0030] After the calculation is completed, the UPDATE module updates the current round and determines whether it is necessary to continue to start the next SMEM search algorithm execution module.

[0031] Preferably, during the SMEM search algorithm calculation process, for a sequence to be aligned, the forward search is calculated sequentially from left to right according to the base order. When it is impossible to continue expanding, the backward search module is started to perform a backward search from right to left starting from the current position. At this time, the forward search module starts a new round of forward search from the point where it cannot be expanded.

[0032] A single loop computation process includes four steps: data initialization, accessing DDR to obtain reference genome data, calculating the backward matching algorithm to update the double interval, and determining whether to continue the loop based on the double interval. Within the search module, the order of computation and memory access is swapped, with computation preceding memory access. Each memory access is for the initialization data required for the next loop, thus achieving a fully pipelined and unblocked loop computation process.

[0033] The double-interval calculation steps in the forward and backward search modules all employ a backward matching algorithm, using a lookup table and an addition tree structure. The four bases A, C, G, and T are encoded with 2 bits. All bases after a specific line in the reference genome are set to A bases. The four bases are grouped together and input into an 8-bit lookup table. The four base codes serve as the lookup table addresses. The output is 32 bits of data representing the number of A, C, G, and T bases in the four bases. The number of each base is represented by 8 bits. Then, the lookup table data results for all bases are accumulated using an addition tree structure. Finally, the number of A bases modified at the beginning is subtracted to obtain the total number of A, C, G, and T bases in the current reference genome substring.

[0034] Preferably, during fractional matrix calculation, the matrix grid points are filled sequentially from the top left to the bottom right diagonal, and the data on the anti-diagonal side are calculated in parallel. The intermediate result data is pulsatingly transmitted in the array through cascaded PE units. In the Smith-Waterman algorithm accelerator, the sequence information storage module adopts a ping-pong cache structure to provide the data required for the current matrix calculation to the pulsating array, while sending the memory access command for the next seed sequence information.

[0035] In the accelerator state machine, the INIT state is responsible for the data initialization of the systolic array. The JUDGE_LEFT state determines whether the current seed can be extended to the left. If it can be extended to the left, it enters the SW state to start the systolic array. Otherwise, it enters the JUDGE_RIGHT state to determine whether it can be extended to the right. If it can be extended to the right, it enters the SW state to start the systolic array again. Otherwise, it enters the FINISH state to complete the bidirectional extension of the current seed.

[0036] In the FINISH state, the score result is written back to DDR and the IDLE state is returned. In the IDLE state, the systolic array hands with the ping-pong buffer of the sequence data storage module. After switching the buffer, the bidirectional extension of the next seed begins. Each time the SW state is entered, the accelerator state machine controller module provides the systolic array with the initialization data and direction of the current extension. The systolic array internally determines the memory read direction and matrix calculation direction based on the initialization data.

[0037] Preferably, the calculation of matrix fractions depends on the input of data in the upper left, upper and left sides, which are represented by H, E and F respectively. During matrix calculation, one PE unit is responsible for the operation of one column of data in the matrix. Therefore, the E fraction in the upper part of the matrix only needs to be maintained in the PE unit, while the other fractions are input through external means. The PE unit has a selector circuit for the input, which selects whether to receive the output result of the previous PE or the initial input according to the control signal. The selection method of other input signals is the same.

[0038] The Smith-Waterman algorithm in BWA-MEM uses a radiometric penalty mechanism. Based on the base similarity between the reference genome bases and the base sequence of the target genome, it updates the H, E, and F scores and outputs them. The PE unit also receives the current maximum score calculated by the previous PE and the matrix column where the maximum score is located, and compares it with the current calculation result. The maximum value is selected and transmitted to the next PE. This score is finally transmitted to the score calculator to calculate the final score of the SW algorithm.

[0039] In the pulsating array, enable signals and initialization data need to be provided periodically and precisely to the PE units, data routers, and fraction calculators. The matrix calculation in the Smith-Waterman algorithm accelerator adopts a bandwidth strategy, with only the diagonal and its surrounding area participating in the calculation. 32 cascaded PEs in the pulsating array complete the matrix filling. In the pulsating array controller, each PE is configured with a PE unit controller, which provides the corresponding PE unit with precise enable signals and determines the matrix column that the PE unit is responsible for calculating.

[0040] A counter is designed in the PE controller circuit. During the counter counting period, the range of enable signals that the controller should send to the PE unit is calculated in advance. The enable signals correspond to the first and last pointers of the count value. When the count value reaches the counter termination point when the controller controls the corresponding PE to calculate the current column, it is determined whether the controlled corresponding PE needs to calculate a new column. If so, the values ​​of the counter start point and termination point, as well as the first and last pointers of the enable range, are updated.

[0041] The heterogeneous acceleration system implementation method for BWA-MEM sequence alignment software provided by the present invention includes the following steps:

[0042] Step 1: The CPU receives the input sequence to be aligned and the compressed file of the reference genome data, parses it to obtain all the sequence data information to be aligned, and then transmits the input data results to the FPGA board to start the SMEM search algorithm accelerator.

[0043] Step 2: Generate seeds for all input sequences in batch processing using the SMEM search algorithm accelerator on the FPGA. After the operation is completed, write the results back to the CPU. The CPU parses the results of the SMEM search algorithm accelerator and then executes the subsequent alignment algorithm for each sequence in sequence.

[0044] Step 3: The CPU program filters all overlapping seeds in a sequence and splices seeds that are close in distance into a seed chain;

[0045] Step 4: The CPU program transmits all the seed chains in a sequence to the FPGA, starts the Smith-Waterman algorithm accelerator on the FPGA, and after the accelerator completes the bidirectional extension of all the seeds in batch processing, the result is written back to the CPU.

[0046] Step 5: Analyze the Smith-Waterman algorithm accelerator results using the CPU to obtain the extension score of each seed, filter them, and save the comparison results.

[0047] Compared with the prior art, the present invention has the following beneficial effects:

[0048] (1) The present invention cuts the BWA-MEM software into fine-grained segments, deploys hot functions to FPGA for acceleration, and the heterogeneous strategy executed by the CPU in the remaining part effectively accelerates sequence alignment. According to the batch processing strategy provided by the present invention, the number of data transmissions between the CPU and FPGA is reduced, the PCIe latency of the heterogeneous system is reduced, and the throughput of the system is improved. It can receive a large number of sequences to be aligned and transmit them to the FPGA for computation in a short time.

[0049] (2) The SMEM search algorithm accelerator provided by this invention decouples the data dependencies between searches, and performs forward and backward searches in parallel pipeline. The search module adopts a fully pipelined architecture, effectively utilizing DDR memory access bandwidth and improving memory access and computation efficiency.

[0050] (3) The Smith-Waterman algorithm accelerator provided by this invention proposes a data reassembly scheme, realizes batch data transmission between CPU and FPGA, uses a pulsating array structure to calculate the fraction matrix, and combines bandwidth strategy and PE reuse strategy. Attached Figure Description

[0051] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0052] Figure 1 A hardware and software block diagram of a CPU and FPGA heterogeneous acceleration system;

[0053] Figure 2 CPU program software reconfiguration diagram;

[0054] Figure 3 This is a diagram of an FPGA-accelerated system architecture.

[0055] Figure 4 Diagram of the SMEM search algorithm accelerator architecture;

[0056] Figure 5 State machine diagram of the SMEM search algorithm;

[0057] Figure 6 This is a pipeline diagram of a base counting circuit;

[0058] Figure 7 This is a diagram of the Smith-Waterman algorithm accelerator architecture.

[0059] Figure 8 State machine diagram of the Smith-Waterman algorithm;

[0060] Figure 9 This is the circuit diagram for the PE unit;

[0061] Figure 10 Diagram of the control strategy for the PE unit controller;

[0062] Figure 11 A graph illustrating the performance of the SMEM search algorithm. Detailed Implementation

[0063] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that those skilled in the art can make several changes and improvements without departing from the concept of the present invention. These all fall within the protection scope of the present invention.

[0064] Example:

[0065] This invention comprises a CPU-FPGA heterogeneous acceleration system for the BWA-MEM gene sequence alignment software. It proposes a heterogeneous segmentation strategy for hotspot functions and a batch processing strategy to reduce PCIe transmission latency due to frequent data transfers between the CPU and FPGA. The invention designs an FPGA acceleration system based on the Intel FPGA acceleration stack platform, mounting an SMEM search algorithm accelerator and a Smith-Waterman algorithm accelerator, and performs fine-grained optimization in the accelerator design to achieve optimal execution speed.

[0066] The heterogeneous acceleration system for BWA-MEM sequence alignment software provided by the present invention includes a CPU program and an FPGA accelerator.

[0067] The CPU program adopts a heterogeneous batch processing strategy to parse the input comparison file, start the FPGA accelerator in sequence, process the FPGA accelerator output, and execute the seed filtering algorithm and seed chain generation algorithm in BWA-MEM.

[0068] The heterogeneous batch processing strategy reconstructs the SMEM search algorithm in the seed generation stage and the Smith-Waterman algorithm in the seed extension stage.

[0069] The SMEM search algorithm executes this part of the program in advance. All input sequences to be compared are transmitted to the FPGA, and the accelerator calculates and obtains the super maximum exact matching seed for all sequences.

[0070] The Smith-Waterman algorithm pre-executes the bidirectional extension Smith-Waterman algorithm, batch-processes all seed chains in a sequence and the bidirectional extension of all seeds within them, and then writes the score results back to the CPU for further filtering.

[0071] Preferably, the FPGA accelerator includes an SMEM search algorithm accelerator;

[0072] The SMEM search algorithm accelerator includes a register configuration and state machine module, a forward search and backward search execution module, and a sequence and intermediate result storage and control module.

[0073] The register configuration and state machine module provides internal initialization data configuration for the accelerator, records the working state of the accelerator, and controls the accelerator's operation process through the state machine.

[0074] The forward search and backward search execution modules decouple the forward and backward search data dependencies within the software, enabling them to run in parallel. Within a single search module, the order of computation and memory access is exchanged, realizing a fully pipelined internal design structure.

[0075] The sequence and intermediate result storage and control module uses BRAM to store the bases of the sequence to be compared and the intermediate results during the search process in a dual-interval manner, and controls the writing and reading of data according to the read and write enable.

[0076] Preferably, the FPGA accelerator also includes a Smith-Waterman algorithm accelerator;

[0077] The Smith-Waterman algorithm accelerator includes a register configuration and state machine module and a systolic array module;

[0078] The pulse array module includes a PE unit array module, a pulse array controller, a data router, and a score calculator module;

[0079] The PE cell array module calculates the score of a matrix grid point within one period according to the Smith-Waterman algorithm score calculation rules.

[0080] The pulse array controller is responsible for providing enable signals to all other modules periodically and precisely.

[0081] The data router is responsible for providing different initialization data to different PE unit array modules, as well as the base data required by the PE unit array modules during operation.

[0082] The score calculator is responsible for recording the output of the score results of the PE unit array module, updating the current maximum score in each cycle, and recording the row and column where the maximum score is located.

[0083] Furthermore, preferred embodiments of the present invention will be described below with reference to the accompanying drawings.

[0084] Figure 1This is a CPU-FPGA heterogeneous system framework. The CPU receives the input alignment sequences and the compressed file of the reference genome data, parses it to obtain all the alignment sequence data, and then transmits the input data results to the FPGA board to start the SMEM search algorithm accelerator. The SMEM search algorithm accelerator on the FPGA batch processes the seed generation of all input sequences. After the operation, the result is written back to the CPU. The CPU parses the SMEM search algorithm accelerator result and then executes the subsequent alignment algorithm for each sequence serially. The CPU program filters all overlapping seeds in a sequence and splices seeds with close distances into seed chains. Then, the CPU program transmits all seed chains in a sequence to the FPGA, starts the Smith-Waterman algorithm accelerator on the FPGA, and after the accelerator batch processes all seeds for bidirectional extension, the result is written back to the CPU. The CPU parses the Smith-Waterman algorithm accelerator result, obtains the extension score of each seed, performs filtering, and saves the alignment results.

[0085] Figure 2 This is the reconstructed CPU software program. The original BWA-MEM software sequentially compares all input sequences in a loop. The SMEM search algorithm is calculated in the `mem_collect_intv()` function, and the bidirectional Smith-Waterman algorithm is calculated in the `ksw_extend2()` function. Combining the proposed batch processing heterogeneous strategy, the two hot functions are executed in advance. The `mem_collect_intv_fpga()` function starts the SMEM search algorithm accelerator on the FPGA to batch process all input sequences to be compared. The `ksw_extend2_fpga()` function starts the Smith-Waterman algorithm accelerator on the FPGA to batch process the bidirectional extension of all seeds in a sequence. This batch processing heterogeneous strategy pre-executes the serially executed parts of the functions within the loop and sends the batch processing to the FPGA accelerator, reducing the number of interactions between the CPU and FPGA during loop execution. Furthermore, batch data transmission has lower latency compared to sequential transmission. In addition, the FPGA accelerator can also improve algorithm acceleration through batch processing sequences.

[0086] Figure 3This is an FPGA acceleration system built on the Intel acceleration stack. The Intel acceleration stack exposes the FPGA acceleration card's interface (such as PCIe and DDR) controller circuitry to a user-defined Acceleration Function (AFU) area through bus encapsulation. It encapsulates the PCIe controller interface as the CCI-P bus protocol and the DDR controller interface as the Avalon-MemoryMap (AVMM) bus protocol. The acceleration system is built based on the corresponding bus protocols and the IP modules provided by the acceleration stack. A DMA module is integrated for data transfer between the CPU and FPGA, and a CCI-P bus 1-to-4 converter bridge is integrated to mount different accelerator modules. The CPU configures the accelerator's internal registers via the CCI-P bus protocol converted from PCIe. The acceleration system mounts an SMEM search algorithm accelerator and a SW algorithm accelerator. To enable the accelerator modules to synthesize higher clock frequencies, the accelerators and peripheral circuit components are clock-isolated using a CCI-P asynchronous bridge, operating at independent clock frequencies. The peripheral circuit components (DMA, converter bridges, etc.) have the same clock frequency as the CCI-P bus protocol. Different modules operate at different times, ensuring that all modules do not access DDR at the same time. Simply add a selector to the interface where the module connects to DDR to select the corresponding functional module at different times.

[0087] Figure 4 This is a design for an SMEM search algorithm accelerator architecture. The CPU configuration of the accelerator's internal registers determines its operating state, including the amount of data in the batch of sequences to be compared, the DDR base address for memory access, and other information. Figure 5 This is the designed accelerator state machine, which controls the seed generation process for all sequences to be aligned. In the SMEM search algorithm accelerator, only a specific number of sequence alignments can be completed at a time. The accelerator pre-determines the number of times the execution module needs to be called to process all current sequences, completing one SMEM search algorithm execution module call from the INIT state to the UPDATE state. In the LOAD state, all sequences to be aligned in the current batch are moved from DDR to the sequence storage module. In the SMEM state, the execution module is started, with the forward and backward search modules starting sequentially, calculating all sequences to be aligned in the current batch in a pipelined manner. After the calculation is complete, the UPDATE module updates the current round and determines whether to continue starting the next SMEM search algorithm execution module.

[0088] like Figure 11In the SMEM search algorithm computation process, the data dependency between the forward and backward searches is decoupled, allowing the forward and backward search modules to start in a pipelined parallel manner. For a sequence to be aligned, the forward search is calculated sequentially from left to right according to the base order. When it cannot be expanded further, the backward search module is started to perform a backward search from right to left starting from the current position. At this time, the forward search module can start a new round of forward search from the point where it cannot be expanded. The single loop computation process includes four steps: data initialization, accessing DDR to obtain reference genome data, calculating the backward matching algorithm to update the double interval, and determining whether to continue the loop based on the double interval. Within the search module, the order of computation and memory access is swapped, with computation preceding memory access. Each memory access retrieves the initialization data required for the next loop, thus achieving a fully pipelined and non-blocking cyclic computation process.

[0089] Figure 6 This is the circuit design for base counting in the reference genome in the backward matching algorithm. The backward matching algorithm is used in both the double-interval calculation steps in the forward and backward search modules. The design uses a lookup table and an addition tree structure. The four bases A, C, G, and T are encoded with 2 bits. All bases after a specific line in the reference genome are set to "A" bases. The four bases are grouped together and input into an 8-bit lookup table. The four base codes are used as the lookup table addresses. The output is 32 bits of data representing the number of A, C, G, and T bases in the four bases. The number of each base is represented by 8 bits. Then, the lookup table data results of all bases are accumulated using an addition tree structure. Finally, the number of "A" bases modified at the beginning is subtracted to obtain the total number of A, C, G, and T bases in the current reference genome substring.

[0090] Figure 7 The designed Smith-Waterman algorithm accelerator architecture also includes a register configuration module and an accelerator state machine control module. Figure 8This is the accelerator state machine. During fractional matrix calculation, the matrix grid points are filled sequentially from the top left to the bottom right diagonal. Data on the anti-diagonal sides can be calculated in parallel, making the systolic array an effective acceleration circuit. Cascaded PE units allow intermediate result data to be systolically transmitted within the array. In the Smith-Waterman algorithm accelerator, the sequence information storage module uses a ping-pong cache structure, providing the systolic array with the data needed for the current matrix calculation while simultaneously sending memory access commands for the next seed sequence. This ping-pong cache structure effectively masks memory access latency through computation. In the accelerator state machine, the INIT state is responsible for initializing the systolic array data. The JUDGE_LEFT state determines whether the current seed can be extended to the left. If it can, it enters the SW state to start the systolic array; otherwise, it enters the JUDGE_RIGHT state to determine whether it can be extended to the right. Similarly, if it can be extended to the right, it enters the SW state again to start the systolic array; otherwise, it enters the FINISH state to complete the bidirectional extension of the current seed. In the FINISH state, the fractional result is written back to DDR, and the system returns to the IDLE state. The systolic array hands over with the ping-pong buffer of the sequence data storage module in the IDLE state. After switching buffers, it begins bidirectional extension of the next seed. Each time it enters the SW state, the accelerator state machine controller module provides the systolic array with the initialization data and direction of the current extension. The systolic array internally determines the memory read direction and matrix calculation direction based on the initialization data.

[0091] Figure 9 This is the circuit design structure of the PE unit in the pulsating array. This module completes the calculation of the data of one grid element in the matrix within one cycle. The calculation of the matrix score depends on the input of the data to its upper left, upper, and left sides, denoted by H, E, and F, respectively. During matrix calculation, one PE unit is responsible for the calculation of one column of data in the matrix. Therefore, the E score of the upper part of the matrix only needs to be maintained within the PE unit, while the remaining scores are input externally. The PE unit has a selector circuit for the input, which selects whether to receive the output result of the previous PE (H_in, F_in) or the initialization input (H_init_in, F_init_in) according to the control signal. The selection method of other input signals is the same, such as inputting reference genome bases (target_init_in). The Smith-Waterman algorithm in BWA-MEM uses a radiometric penalty mechanism to update the H, E, and F scores and output them based on the base similarity between the reference genome bases and the target genome sequence. In addition, the PE unit will receive the current maximum score calculated by the previous PE and the matrix column where the maximum score is located, compare it with the current calculation result, select the maximum value and transmit it to the next PE. This score is finally transmitted to the score calculator to calculate the final score of the SW algorithm.

[0092] In the pulsating array, enable signals and initialization data need to be periodically and precisely provided to the PE units, data routers, and fraction calculators. The matrix computation in the Smith-Waterman algorithm accelerator adopts a bandwidth strategy, with only the diagonal and its vicinity participating in the computation to reduce the amount of matrix calculation. The matrix width and height of each seed Smith-Waterman algorithm are uncertain. To make the accelerator more generalizable, 32 PEs are cascaded in the pulsating array to complete the matrix filling. In the pulsating array controller, each PE is configured with a PE unit controller, which precisely provides enable signals to the corresponding PE unit and determines the matrix column that the PE unit is responsible for calculating.

[0093] Figure 10 This describes the dual-pointer control strategy of the PE unit controller, using an example of a bandwidth setting of 6 and cascading 8 PE units. The gray shaded areas in the matrix represent the grid points that need to be calculated. Each PE unit is responsible for the calculation of one column of data in the matrix, and the score for each grid point is calculated within one cycle. A counter is designed in the PE controller circuit. Figure 10 In the diagram, S and E represent the start and end points of the counter when the controller controls the corresponding PE to calculate the current column. During the counter counting period, the range of enable signals that the controller should send to the PE unit is calculated in advance, which corresponds to the first and last pointers of the count value. When the count value reaches E, it is determined whether the controlled PE needs to calculate a new column. If so, the values ​​of S and E, as well as the first and last pointers of the enable range, are updated. Figure 10 In this configuration, PE0 is initially responsible for the calculation of the first column. The counter S equals 1, and E equals 11. When the count value is between 1 and 6, the PE controller should send an enable signal to the PE unit. Once the counter value reaches 11, PE0 will be responsible for calculating the elements of the ninth column of the matrix, at which point S and E are updated, along with the enabled count value range. The control methods for the remaining PEs are the same; only different start and end points need to be configured.

[0094] This invention provides a heterogeneous acceleration system implementation method for BWA-MEM sequence alignment software, comprising the following steps:

[0095] Step 1: The CPU receives the input sequence to be aligned and the compressed file of the reference genome data, parses it to obtain all the sequence data information to be aligned, and then transmits the input data results to the FPGA board to start the SMEM search algorithm accelerator.

[0096] Step 2: Generate seeds for all input sequences in batch processing using the SMEM search algorithm accelerator on the FPGA. After the operation is completed, write the results back to the CPU. The CPU parses the results of the SMEM search algorithm accelerator and then executes the subsequent alignment algorithm for each sequence in sequence.

[0097] Step 3: The CPU program filters all overlapping seeds in a sequence and splices seeds that are close in distance into a seed chain;

[0098] Step 4: The CPU program transmits all the seed chains in a sequence to the FPGA, starts the Smith-Waterman algorithm accelerator on the FPGA, and after the accelerator completes the bidirectional extension of all the seeds in batch processing, the result is written back to the CPU.

[0099] Step 5: Analyze the Smith-Waterman algorithm accelerator results using the CPU to obtain the extension score of each seed, filter them, and save the comparison results.

[0100] Those skilled in the art will understand that, in addition to implementing the system, apparatus, and their modules provided by this invention in purely computer-readable program code, the same program can be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, the system, apparatus, and their modules provided by this invention can be considered a hardware component, and the modules included therein for implementing various programs can also be considered structures within the hardware component; alternatively, modules for implementing various functions can be considered both software programs implementing the method and structures within the hardware component.

[0101] Specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Unless otherwise specified, the embodiments and features described in this application can be arbitrarily combined with each other.

Claims

1. A heterogeneous acceleration system for BWA-MEM sequence alignment software, characterized in that, Including CPU programs and FPGA accelerators; The CPU program adopts a heterogeneous batch processing strategy to parse the input comparison file, start the FPGA accelerator in sequence, process the FPGA accelerator output, and execute the seed filtering algorithm and seed chain generation algorithm in BWA-MEM. The heterogeneous batch processing strategy reconstructs the SMEM search algorithm in the seed generation stage and the Smith-Waterman algorithm in the seed extension stage. The SMEM search algorithm executes this part of the program in advance. All input sequences to be compared are transmitted to the FPGA, and the accelerator calculates and obtains the super maximum exact matching seed for all sequences. The Smith-Waterman algorithm is described above. The bidirectional extension Smith-Waterman algorithm is executed in advance. After batch processing all seed chains and all seeds within a sequence, the bidirectional extension is performed, and the score results are written back to the CPU for further filtering. The FPGA accelerator includes an SMEM search algorithm accelerator; The SMEM search algorithm accelerator includes a register configuration and state machine module, a forward search and backward search execution module, and a sequence and intermediate result storage and control module. The register configuration and state machine module provides internal initialization data configuration for the accelerator, records the working state of the accelerator, and controls the accelerator's operation process through the state machine. The forward search and backward search execution modules decouple the forward and backward search data dependencies within the software, enabling them to run in parallel. Within a single search module, the order of computation and memory access is swapped to achieve a fully pipelined internal design structure. The sequence and intermediate result storage and control module uses BRAM to store the bases of the sequence to be compared and the intermediate results during the search process in a dual-interval manner, and controls the writing and reading of data according to the read and write enable. The FPGA accelerator also includes a Smith-Waterman algorithm accelerator; The Smith-Waterman algorithm accelerator includes a register configuration and state machine module and a systolic array module; The pulse array module includes a PE unit array module, a pulse array controller, a data router, and a score calculator module; The PE cell array module calculates the score of a matrix grid point within one period according to the Smith-Waterman algorithm score calculation rules. The pulse array controller is responsible for providing enable signals to all other modules periodically and precisely. The data router is responsible for providing different initialization data to different PE unit array modules, as well as the base data required by the PE unit array modules during operation. The score calculator is responsible for recording the output of the score results of the PE unit array module, updating the current maximum score in each cycle, and recording the row and column where the maximum score is located.

2. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, wherein, Two hot functions are executed in advance: in the mem_collect_intv_fpga() function, the SMEM search algorithm accelerator on the FPGA is started to batch process all input sequences to be compared; in the ksw_extend2_fpga() function, the Smith-Waterman algorithm accelerator on the FPGA is started to batch process the bidirectional extension of all seeds in a sequence.

3. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, characterized in that, The Intel acceleration stack is used to encapsulate the interface controller circuit of the FPGA acceleration card into a bus and expose it to the user-defined acceleration function area. The PCIe controller interface is encapsulated into the CCI-P bus protocol, and the DDR controller interface is encapsulated into the AVMM bus protocol. An acceleration system is built based on the corresponding bus protocol and the IP modules provided by the acceleration stack. A DMA module is integrated for data transmission between the CPU and FPGA. A CCI-P bus one-to-four-converter bridge is integrated for mounting different accelerator modules. The CPU configures the internal registers of the accelerator through the CCI-P bus protocol converted from PCIe. The acceleration system is equipped with an SMEM search algorithm accelerator and a SW algorithm accelerator. The accelerators and peripheral circuit components are clock isolated by a CCI-P asynchronous bridge and operate at an independent clock frequency. The peripheral circuit components have the same clock as the CCI-P bus protocol.

4. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, characterized in that, In the SMEM search algorithm accelerator, only a specific number of sequence alignments are completed at a time. The accelerator obtains in advance the number of times the execution module needs to be called in a loop to execute all current sequences, and completes one SMEM search algorithm execution module call from the INIT state to the UPDATE state. In the LOAD state, all the sequences to be compared that need to be processed in the current batch are moved from DDR to the sequence storage module. In the SMEM state, the execution module is started and the forward and backward search modules are started in sequence to calculate all the sequences to be compared in the current batch in a pipeline manner. After the calculation is completed, the UPDATE module updates the current round and determines whether it is necessary to continue to start the next SMEM search algorithm execution module.

5. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, characterized in that, During the SMEM search algorithm calculation process, for a sequence to be aligned, the forward search is calculated sequentially from left to right according to the base order. When it is impossible to continue expanding, the backward search module is started to perform a backward search from right to left starting from the current position. At this time, the forward search module starts a new round of forward search from the point that cannot be expanded. A single loop computation process includes four steps: data initialization, accessing DDR to obtain reference genome data, calculating the backward matching algorithm to update the double interval, and determining whether to continue the loop based on the double interval. Within the search module, the order of computation and memory access is swapped, with computation preceding memory access. Each memory access is for the initialization data required for the next loop, thus achieving a fully pipelined and unblocked loop computation process. The double-interval calculation steps in the forward and backward search modules all employ a backward matching algorithm, using a lookup table and an addition tree structure. The four bases A, C, G, and T are encoded with 2 bits. All bases after a specific line in the reference genome are set to A bases. The four bases are grouped together and input into an 8-bit lookup table. The four base codes serve as the lookup table addresses. The output is 32 bits of data representing the number of A, C, G, and T bases in the four bases. The number of each base is represented by 8 bits. Then, the lookup table data results for all bases are accumulated using an addition tree structure. Finally, the number of A bases modified at the beginning is subtracted to obtain the total number of A, C, G, and T bases in the current reference genome substring.

6. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, characterized in that, When calculating the fractional matrix, the matrix grid points are filled sequentially from the top left to the bottom right diagonal, and the data on the anti-diagonal are calculated in parallel. The intermediate result data is pulsatingly transmitted in the array through cascaded PE units. In the Smith-Waterman algorithm accelerator, the sequence information storage module adopts a ping-pong cache structure to provide the data needed for the current matrix calculation to the pulsating array, and at the same time send the memory access command for the next seed sequence information. In the accelerator state machine, the INIT state is responsible for the data initialization of the systolic array. The JUDGE_LEFT state determines whether the current seed can be extended to the left. If it can be extended to the left, it enters the SW state to start the systolic array. Otherwise, it enters the JUDGE_RIGHT state to determine whether it can be extended to the right. If it can be extended to the right, it enters the SW state to start the systolic array again. Otherwise, it enters the FINISH state to complete the bidirectional extension of the current seed. In the FINISH state, the score result is written back to DDR and the IDLE state is returned. In the IDLE state, the systolic array hands with the ping-pong buffer of the sequence data storage module. After switching the buffer, the bidirectional extension of the next seed begins. Each time the SW state is entered, the accelerator state machine controller module provides the systolic array with the initialization data and direction of the current extension. The systolic array internally determines the memory read direction and matrix calculation direction based on the initialization data.

7. The heterogeneous acceleration system for BWA-MEM sequence alignment software according to claim 1, characterized in that, The calculation of matrix fractions depends on the input of its upper left, upper, and left sides, respectively, using... , and This means that during matrix computation, one PE unit is responsible for the operation on one column of data in the matrix; therefore, the area above the matrix... The fractions only need to be maintained within the PE unit, while the remaining fractions are obtained through external input. The PE unit has a selector circuit for the input, which selects whether to receive the output result of the previous PE or the initial input based on the control signal. The selection method for the remaining input signals is the same. The Smith-Waterman algorithm in BWA-MEM employs a radiation penalty mechanism, updating the base similarity between the reference genome and the target genome sequence based on the base similarity. , and The PE unit also receives the current maximum score calculated by the previous PE and the matrix column where the maximum score is located, compares it with the current calculation result, selects the maximum value and transmits it to the next PE. This score is finally transmitted to the score calculator to calculate the final score of the SW algorithm. In the pulsating array, enable signals and initialization data need to be provided periodically and precisely to the PE units, data routers, and fraction calculators. The matrix calculation in the Smith-Waterman algorithm accelerator adopts a bandwidth strategy, with only the diagonal and its surrounding area participating in the calculation. 32 cascaded PEs in the pulsating array complete the matrix filling. In the pulsating array controller, each PE is configured with a PE unit controller, which provides the corresponding PE unit with precise enable signals and determines the matrix column that the PE unit is responsible for calculating. A counter is designed in the PE controller circuit. During the counter counting period, the range of enable signals that the controller should send to the PE unit is calculated in advance. The enable signals correspond to the first and last pointers of the count value. When the count value reaches the counter termination point when the controller controls the corresponding PE to calculate the current column, it is determined whether the controlled corresponding PE needs to calculate a new column. If so, the values ​​of the counter start point and termination point, as well as the first and last pointers of the enable range, are updated.

8. A method for implementing a heterogeneous acceleration system for BWA-MEM sequence alignment software, characterized in that, The heterogeneous acceleration system for BWA-MEM sequence alignment software according to any one of claims 1-7 includes the following steps: Step 1: The CPU receives the input sequence to be aligned and the compressed file of the reference genome data, parses it to obtain all the sequence data information to be aligned, and then transmits the input data results to the FPGA board to start the SMEM search algorithm accelerator. Step 2: Generate seeds for all input sequences in batch processing using the SMEM search algorithm accelerator on the FPGA. After the operation is completed, write the results back to the CPU. The CPU parses the results of the SMEM search algorithm accelerator and then executes the subsequent alignment algorithm for each sequence in sequence. Step 3: The CPU program filters all overlapping seeds in a sequence and splices seeds that are close in distance into a seed chain; Step 4: The CPU program transmits all the seed chains in a sequence to the FPGA, starts the Smith-Waterman algorithm accelerator on the FPGA, and after the accelerator completes the bidirectional extension of all the seeds in batch processing, the result is written back to the CPU. Step 5: Analyze the Smith-Waterman algorithm accelerator results using the CPU to obtain the extension score of each seed, filter them, and save the comparison results.

Images