Method for gene sequence alignment and electronic equipment

Abstract

The invention provides a gene sequence alignment method and electronic equipment, the alignment process is divided into a seed search stage and a seed expansion stage, in the seed search stage, a plurality of seeds are determined for a read sequence to be aligned; acquiring a plurality of query sequences to be compared and corresponding reference sequences thereof based on the seeds in a seed expansion stage; expanding the length of each sequence by an integral multiple of the data width of the SIMD instruction, and initializing a score matrix; in the calculation process of each round of comparison score, the calculation range of the current round is determined in each score matrix, and then the comparison score is calculated by using an SIMD comparison instruction in the reverse diagonal direction of the score matrix; and after all rounds of calculation are completed, determining an optimal matching sequence according to the maximum comparison score. According to the scheme, the sequence length difference in the comparison process is effectively considered, a large amount of invalid calculation is reduced, the operation speed of the seed expansion stage is increased, and then the sequence comparison efficiency is improved.

Description

Methods and electronic devices for gene sequence alignment Technical Field

[0001] This application relates to biological sequencing and gene data analysis, and in particular to methods and apparatus for gene sequence alignment. Background Technology

[0002] With the development of precision medicine and gene sequencing technology, the amount of gene data used for computation is growing exponentially, and the efficiency of gene data analysis is receiving increasing attention. Whether it's second-generation short-read sequencing data or third-generation long-read sequencing data, gene sequence alignment is a very time-consuming step. Therefore, how to improve the efficiency of gene sequence alignment has become a research hotspot.

[0003] Current mainstream gene sequence alignment technologies generally follow a "seed-expansion" model. First, in the seed search stage, precise matching techniques quickly identify short fragments (called "seeds") that are highly similar and perfectly matched between the read sequence to be aligned (also called the query sequence) and the reference genome (also called the reference sequence). Then, in the seed expansion stage, the sequence extends and aligns base by base from the location of these seeds (allowing mismatches, insertions / deletions), calculating alignment scores and selecting the alignment with the highest score as the output (i.e., the best match). The seed expansion stage requires calculating the expansion score of each seed on the reference genome, which is computationally intensive. Improving the performance of the seed expansion stage in the alignment process has become a key challenge in gene sequence alignment. Summary of the Invention

[0004] This application aims to improve the speed of the seed expansion stage in gene sequence alignment, thereby improving the efficiency of gene sequence alignment.

[0005] The objective of this application is achieved through the following technical solution:

[0006] According to a first aspect of the embodiments of this application, a method for gene sequence alignment is provided, comprising a seed search stage and a seed expansion stage, wherein the seed expansion stage comprises: obtaining multiple query sequences to be aligned and their corresponding reference sequences based on multiple seeds determined by the seed search stage for read sequences to be aligned; expanding the length of these sequences to an integer multiple of the data width of the SIMD instruction and initializing a score matrix for each query sequence and its corresponding reference sequence; determining the calculation range of the current round in the score matrix, and calculating the alignment score of the query sequence and its corresponding reference sequence using the SIMD comparison instruction along the anti-diagonal direction of the score matrix; in response to the maximum alignment score of the current round being greater than the maximum alignment score of the previous round, obtaining the alignment position corresponding to the maximum alignment score of the current round; and after all rounds have been calculated, determining the best matching sequence between the query sequence and the reference sequence based on the alignment position corresponding to the maximum alignment score and the score matrix.

[0007] In this embodiment, the query and reference sequences participating in the alignment are padded according to the data width of the SIMD instructions. In each round of iteration, the alignment score is calculated along the anti-diagonal of the score matrix. This eliminates the need to wait for the last sequence to be calculated before proceeding to the next round, thus reducing idle and wasted computational resources when shorter query sequences are compared to longer ones. Furthermore, if the maximum alignment score in the current round does not exceed the maximum alignment score in the previous round, the calculation of the coordinates corresponding to the maximum alignment score in the current round is ignored, saving the need for iterative calculations. Therefore, this embodiment improves the efficiency of gene sequence alignment, especially suitable for situations with unbalanced query sequence lengths.

[0008] In some embodiments, the method may further include: determining the maximum alignment score for the current round after the alignment scores for the current round have been calculated by: performing an element-wise maximum maximization operation on the alignment scores already calculated for the current round, arranged along the anti-diagonal direction of the score matrix, and recording each maximum value in the sliding window data structure, wherein the length and step size of the sliding window data structure are equal to the data width of the SIMD instruction; and moving the sliding window data structure along the anti-diagonal direction to traverse all alignment scores for the current round, thereby determining the maximum alignment score for the current round. In this embodiment, the length of the sliding window data structure is set to the data width of the SIMD instruction to facilitate the execution of the SIMD instruction. Using the SIMD MAX instruction in parallel to determine the maximum alignment score for each round reduces the number of instructions used, thereby improving computational efficiency.

[0009] In some embodiments, the direction of moving the sliding window data structure is from the lower left corner to the upper right corner.

[0010] In some embodiments, determining the calculation range for the current round may include determining the calculation range based on a pre-set calculation width parameter, which is set based on the matching position of the query sequence and the reference sequence. Before each round of calculation, the range that needs to be calculated within the score matrix is ​​reduced according to the calculation width parameter, instead of directly calculating every element on the diagonal line of that round. This reduces the amount of computation and avoids affecting high-value areas in the score matrix, thereby ensuring the correctness of the final obtained best matching sequence.

[0011] In some embodiments, the method may further include: after the calculation of the current round is completed, traversing the comparison scores of the current round and the previous round along both ends of the anti-diagonal direction; if, at either end, the comparison scores of the current round and the previous round simultaneously show consecutive zero values, then the calculation range is narrowed down to the non-zero score positions; the narrowed range is compared with the calculation width parameter, and the smaller one is selected as the calculation range for the next round. In this embodiment, the above-mentioned method of further narrowing the calculation range (or pruning) identifies the low-value regions in the comparison scores to be calculated in the next round. By pruning these regions, the number of actual comparison scores calculated is reduced, thus improving calculation efficiency while ensuring the comparison results.

[0012] In some embodiments, the method may further include: when the calculation range narrows to 0, stopping the calculation of the next round and subsequent rounds, and using the alignment position corresponding to the maximum alignment score of the current record as the score matrix to determine the best matching sequence between the query sequence and the reference sequence. When the above situation occurs, it indicates that the next round and subsequent rounds are both low-value areas and do not require calculation. By stopping the calculation of subsequent rounds in a timely manner through this judgment, invalid calculations can be further avoided.

[0013] In some embodiments, after determining the calculation range of the current round in the score matrix, the query sequence and reference sequence containing the parts to be compared are read using SIMD read instructions. Then, the data width of the SIMD instructions is processed simultaneously in each calculation, thereby realizing parallel calculation of the comparison scores.

[0014] In some embodiments, calculating the alignment score between the query sequence and its corresponding reference sequence includes: when the calculated element is within the calculation range of this round in the score matrix, calculating using the following formula:

[0015]

[0016] in, For the reference sequence The element and the query sequence of the first element The matching score of each element. The penalty points are calculated based on the penalty score.

[0017] In some embodiments, determining the best matching sequence between the query sequence and the reference sequence includes backtracking from the maximum alignment score in the score matrix along the upper left, left, or upper direction to the element with an alignment score of 0. In this embodiment, the best matching sequence is determined by this backtracking method, and the relationship between the bases of the two gene sequences is also determined simultaneously, such as matching (or substitution), insertion, or deletion.

[0018] According to a second aspect of the present application, an electronic device is provided, the electronic device including a memory and a processor, the memory storing computer-executable instructions executed by the processor, the computer-executable instructions including SIMD instructions, the processor executing the computer-executable instructions to implement the method described in the first aspect of the present application.

[0019] According to a third aspect of the present application, a computer-readable storage medium is provided that stores a computer program or computer instructions, which are executed by a processor to implement the method described in the first aspect of the present application.

[0020] According to a fourth aspect of the embodiments of this application, a computer program product is provided, which includes computer program code or computer instructions, and when the computer program code or computer instructions are run, implements the method described in the first aspect of the embodiments of this application. Attached Figure Description

[0021] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application. It is obvious that the drawings described below are merely some embodiments of this application, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort. In the drawings:

[0022] Figure 1 shows a schematic flowchart of a gene sequence alignment method according to an embodiment of this application;

[0023] Figure 2 shows a schematic diagram of each round of calculation of the comparison score in the example;

[0024] Figure 3 shows a schematic diagram of each round of calculation of the comparison score in the example;

[0025] Figure 4 illustrates an example of the process for determining the maximum comparison score in each round; and

[0026] Figure 5 shows a schematic diagram of a dynamic pruning process according to an embodiment of this application. Detailed Implementation

[0027] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description, in conjunction with the accompanying drawings and specific embodiments, further illustrates this application. It should be understood that the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application.

[0028] Furthermore, the described features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. Numerous specific details are provided in the following description to give a thorough understanding of embodiments of this application. However, those skilled in the art will recognize that the technical solutions of this application can be practiced without one or more of the specific details, or other methods, components, apparatuses, steps, etc., can be employed. In other instances, well-known methods, apparatuses, implementations, or operations are not shown or described in detail to avoid obscuring various aspects of this application.

[0029] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0030] The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be broken down, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0031] The mainstream sequence alignment method currently employs a seed-extension approach, comprising two stages: a seed lookup stage and a seed extension stage. In the seed lookup stage, multiple short sequence fragments (seeds) that perfectly match the query sequence and the reference sequence are quickly identified. In the seed extension stage, bidirectional extension is performed from the seed location (also called the seed site) along the start and end coordinates of the seed matching, extracting corresponding fragments from the query and reference sequences to form local sequence pairs to be aligned. Then, a local alignment algorithm (such as the Smith-Waterman algorithm) is used to align the extended local sequence pairs, calculate the alignment score, and select the optimal alignment result.

[0032] The Smith-Waterman algorithm (SW algorithm), proposed by T.S. Smith and MS Waterman in 1981, is a local alignment algorithm based on dynamic programming. Its goal is to find highly similar sequence segments between two sequences. The main steps include: constructing a score matrix, backtracking the optimal path, and generating the alignment results. The SW algorithm primarily relies on the alignment score matrix to determine the maximum alignment score. For example, each element in the alignment score matrix of sequence A and sequence B... Represents the preceding sequence A The characters and sequence B before The optimal local alignment score for each character is determined by backtracking based on the matrix of alignment scores and the maximum alignment score, thereby determining the optimal matching sequence.

[0033] The formula for calculating each element in the basic score matrix is ​​as follows:

[0034] (1)

[0035] in, For the reference sequence elements (i.e.) ) and the query sequence elements (i.e.) The matching score of ) As a penalty score, The length of the reference sequence is greater than or equal to 1. The length of the reference sequence is greater than or equal to 1 and less than or equal to the length of the query sequence. According to this formula, when the length of the reference sequence is known... The element and the query sequence of the first element When there are 1 element, it can be based on the matrix elements. , and The value of the matrix element is used to calculate the matrix element. The value of .

[0036] In existing sequence alignment software, each seed is not typically expanded independently. Instead, preprocessing such as clustering and merging is performed on the seeds to determine multiple candidate matching regions. Then, a fixed length is extended outwards from the candidate matching region to define the local window for Smith-Waterman (SW) alignment (i.e., the reference sequence window and its corresponding query sequence window). Based on this, a score matrix is ​​constructed, with the dimension of the query sequence window length × the reference sequence window length. The SW score matrix is ​​forced to calculate only the region within ±k bp of the diagonal (k is usually 10~20 bp). This approach is also known as the Banded Smith-Waterman (BSW) algorithm, which limits the actual calculation range to a single band within the entire matrix, allowing for rapid determination of the maximum alignment score.

[0037] In improved sequence alignment software (such as BWA-MEM2), an inter-query parallel computing method is further introduced on the basis of BSW. By using SIMD instructions, the computational efficiency of the seed expansion stage is improved, and multiple sequences can be processed simultaneously.

[0038] However, through in-depth research and practice, the inventors discovered that this technique of accelerating gene alignment using SIMD instructions still has many limitations, especially in scenarios where the lengths of multiple reference and query sequences involved in the alignment vary. For ease of description, the two sequences involved in the SW local alignment will be referred to as the reference sequence and the query sequence in the following text. BWA-MEM2 assigns multiple independent query sequences to independent vector slots in the SIMD register. Without changing the row-by-row left-to-right calculation direction and logic of the SW algorithm for a single sequence, a single SIMD instruction performs the same SW calculation step (which can be understood as one round) simultaneously on all sequences in all vector slots, thereby achieving parallel acceleration of multiple sequences. For example, a sequence of length... The matrix corresponding to the query sequence is Another length is The matrix corresponding to the query sequence is In the The wheels need to be counted from left to right. … as well as … ,calculate … With calculation … It is parallel. However, according to formula (1) above, the calculation... need The value, calculate need The value of is calculated, and so on, therefore in the calculation … (or … When performing alignment, calculations must be performed sequentially and cannot be parallelized. In other words, this method requires multiple rounds of parallel computation on these sequences to complete the entire SW alignment process. However, for each sequence, the alignment score calculated in each round (except for the first calculation in the current round) requires the result of the previous calculation in the current round. For example, calculating... need The value, therefore, the result that needs to be calculated in the same round, such as and It cannot be computed in parallel.

[0039] The inventors discovered that in this comparison process, if the lengths of the various query sequences differ significantly, in each round of calculation, it is necessary to wait for each query sequence to complete its calculation before proceeding to the next round. However, by this time, the shorter query sequence has already finished its calculation, forcing the inventors to wait for the longer query sequence to complete its round. For example, in the query sequences mentioned above, if... Greater than Then the length is The query sequence completes this round of calculation and calculates to After that, the length is The query sequence may only be calculated up to Not yet calculated However, this round of calculations was not completed, therefore it was used for a length of The computational resources for the query sequence will be in the length of The query sequence remains idle while computation continues. This wastes computational resources, thus reducing the efficiency of gene sequence alignment.

[0040] For scenarios where the lengths of the reference and query sequences vary, this application provides a novel sequence alignment scheme that can more effectively utilize SIMD vector instructions, reduce a large amount of invalid computation, improve computational efficiency, accelerate the operation of the seed expansion stage, and thus improve the efficiency of sequence alignment.

[0041] Figure 1 shows a schematic flowchart of a gene sequence alignment method 100 according to an embodiment of this application. The method 100 includes: S110, obtaining multiple query sequences to be aligned and their corresponding reference sequences based on multiple seeds determined for the read sequences to be aligned via a seed search phase; S120, expanding the length of these sequences to an integer multiple of the data width of the SIMD instruction and initializing a score matrix for each query sequence and its corresponding reference sequence; S130, determining the calculation range of the current round in the score matrix, and calculating the alignment score of the query sequence and its corresponding reference sequence using the SIMD comparison instruction along the anti-diagonal direction of the score matrix; S140, in response to the maximum alignment score of the current round being greater than the maximum alignment score of the previous round, obtaining the alignment position corresponding to the maximum alignment score of the current round; S150, after all rounds of calculation are completed, determining the best matching sequence between the query sequence and the reference sequence based on the alignment position corresponding to the maximum alignment score and the score matrix.

[0042] More specifically, in S110, centered on the location of each seed determined in the seed search phase, a bidirectional expansion is performed from the start and end coordinates of seed matching, extracting corresponding fragments from the read sequences and the reference genome, thereby forming local sequence pairs to be aligned, i.e., multiple query sequences and their corresponding reference sequences. In subsequent processing, an alignment score is calculated for each query sequence and its corresponding reference sequence.

[0043] In S120, the lengths of the acquired multiple query sequences and their corresponding reference sequences are first expanded to an integer multiple of the data width of the SIMD instruction, with any portion exceeding the original sequence assigned a zero value. This not only aligns the data, facilitating the allocation of contiguous storage space for data storage, but also fully utilizes Single Instruction Multiple Data (SIMD) vectorization technology. After expansion, the lengths of the data read, computed, and written by the SIMD instruction are the same, and computations are performed on vectors of the same dimension, increasing the parallelism of SIMD instruction execution and improving computational efficiency.

[0044] Then, a score matrix is ​​initialized for each query sequence and its corresponding reference sequence. Each score matrix is ​​a two-dimensional matrix, with the number of rows and columns related to the length of the query sequence and the length of its corresponding reference sequence. For example, if the length of the expanded reference sequence is m and the length of the expanded query sequence is n, then the corresponding score matrix has m+1 rows and n+1 columns. The 0th row and 0th column of the score matrix do not represent the highest local alignment score between a base in the reference sequence and a base in the query sequence, but are used for the calculation of the 1st row and 1st column of the score matrix. The other elements of the score matrix, such as... , indicating the preceding part of the reference sequence The first element and query sequence The highest local alignment score of each element.

[0045] In S130, SIMD instructions are used to simultaneously compute multiple acquired query sequences and their reference sequences. As mentioned above, the local alignment algorithm SW is completed through calculations in each round. In each round, the same SIMD instruction is used to perform the same processing on all sequences simultaneously, thereby achieving parallel acceleration of multiple sequences. Before the calculation begins in each round, the calculation range for the current round is first determined in the score matrix. For example, the calculation range can be determined according to a pre-set calculation width parameter, which is set based on the matching position of the query sequence and the reference sequence. A smaller width results in faster computation, but the approximate maximum alignment score calculated may deviate significantly from the actual maximum alignment score, potentially missing the best matching sequence. A larger width, on the other hand, is more sensitive to computation and more likely to obtain an approximate maximum alignment score close to the actual maximum alignment score, but the computational cost will increase significantly. The width parameter can be pre-set according to actual needs. Next, SIMD read instructions are used to read the data containing the query sequence and reference sequence to be compared. These sequence data with the same width can be stored in contiguous storage space.

[0046] During the calculation process, the calculation direction in each round is along the anti-diagonal direction of the score matrix. The anti-diagonal direction includes the direction from the lower left corner to the upper right corner and the direction from the upper right corner to the lower left corner. In a preferred embodiment, the calculation direction in each round is along the direction from the lower left corner to the upper right corner of the score matrix, that is, in each round according to... , , … calculate in the following order (with the vertical downward direction as the x-axis and the horizontal rightward direction as the y-axis).

[0047] Figure 2 illustrates the calculation of the comparison score in each round. Rounds N-1 and N have been completed, and round N+1 is awaiting calculation. As shown, according to the Smith-Waltman algorithm, when the calculation direction in each round is along the anti-diagonal of the score matrix, each element in round N+1 can be calculated based on the results of rounds N-1 and N, and each element in round N+1 can be calculated in parallel without depending on any other element in round N+1. For example, position A in the figure... It can be calculated based on the (N-1)th round. and the result calculated in the Nth round and The calculation is performed without relying on any other elements calculated in the (N+1)th round. (The element at position B in the diagram is shown.) It can be calculated based on the (N-1)th round. The calculation in the Nth round And outside the side strip The alignment score (treated as 0 or assigned a value of 0) is calculated without relying on any other elements calculated in round N+1. Therefore, the alignment scores calculated in each round (e.g., all elements calculated in round N+1 here) can be computed in parallel. It can be seen that although multiple rounds of parallel computation are still required to complete the entire SW alignment process for these sequences, for each sequence, the alignment score calculated in each round no longer requires the result of the previous round, thus avoiding mutual waiting between sequences of varying lengths within the same round.

[0048] For example, Figure 3 shows a schematic diagram of each round of comparison score calculation. Assume the reference sequence R (“AATCGC”) is arranged vertically, the query sequence Q (“ACAGTA”) is arranged horizontally, and the calculation direction for each round is the anti-diagonal direction from the lower left to the upper right. First, the values ​​in row 0 and column 0 of the first region are set, as shown in Figure 1, for subsequent comparison score calculation. In some embodiments, the preset values ​​are all “0”. In this embodiment, the preset values ​​are “7, 1, 0, 0”. The preset values ​​can be determined according to the specific algorithm design and are not limited here. Then, the first round of calculation is performed. At this time, only the following calculation is needed: The value of is compared between A in the reference sequence R and A in the query sequence Q. Since they are the same at this point, the Smith-Waltman algorithm is used to calculate . ,in It is 7. The table shows 1, the result is 8, while other calculation methods all result in 0. Therefore... The value is 8, as shown in the second image. Then, the second round of calculations is performed. At this point, two of the four elements along the diagonal line need to be calculated, namely... and According to the Smith-Waltman algorithm, the matrix shown in Figure 3 is obtained. The third round of calculation is then performed. At this point, the diagonal line is among the three elements in the first region, i.e. , and All need to be calculated, and and Since it is not within the first region, it is not necessary to actually perform calculations according to the Smith-Waltman algorithm. The calculation results are shown in Figure 4.

[0049] In some embodiments, after determining the calculation range for the current round in the scoring matrix, the data width of the SIMD instructions can be processed in parallel in one go. If the calculated data exceeds the data to be compared, the result for the excess data is assigned a value of 0.

[0050] For example, in a certain round, it is necessary to calculate the score matrix. , , , and The dataset has five elements. The SIMD instruction has a data width of 4, and the computation direction is from the bottom left to the top right. In this round of computation, because SIMD instructions are used, it actually involves two computations; for example, the first computation can be performed in parallel. , , and The second one can be computed in parallel. , , and .in, , , , and These five elements are assigned values ​​based on the results of the Smith-Waltman algorithm, and , and For data that does not require comparison, the result is directly assigned the value 0.

[0051] Referring back to Figure 1, in S140, after each round of calculation, the alignment position of the maximum alignment score of the calculated portion needs to be updated. If the maximum alignment score of the current round is greater than the maximum alignment score of the previous round, the alignment position corresponding to the maximum alignment score of the calculated portion is updated, i.e., the alignment position corresponding to the maximum alignment score of the current round is obtained as the alignment position of the maximum alignment score of the calculated portion. If the maximum alignment score of the current round is less than the maximum alignment score of the previous round, it means that the maximum alignment score of a certain previous round is the maximum alignment score of the calculated portion, so there is no need to update the alignment position of the maximum alignment score of the calculated portion, thus avoiding unnecessary traversal. In S150, when all rounds are completed, the maximum alignment score and its alignment position (the corresponding positions of the reference sequence and the query sequence) of the entire score matrix can be obtained. Then, backtracking is performed based on the matrix of calculated alignment scores and the maximum alignment score to determine the optimal matching sequence for the query sequence.

[0052] Figure 4 illustrates a schematic diagram of the process for determining the maximum comparison score in each round according to an embodiment of this application. The figure shows the calculation process for the maximum comparison score in rounds N and N+2:

[0053] First, the calculated alignment scores for this round are subjected to an element-wise max operation using a sliding window data structure. For round N, seven alignment scores have been calculated within the computational range: "0, 5, 6, 9, 6, 0, 4". The length of the sliding window data structure is equal to the length of four alignment scores. In step 0, the sliding window data structure is initialized to "0, 0, 0, 0". In step 1, the sliding window is placed at the initial position of the calculated alignment scores for this round (corresponding to "0, 5, 6, 9"). Then, each element of the sliding window data structure is compared with the corresponding element of the calculated alignment scores for this round. If an element of the sliding window data structure is smaller than an element of the calculated alignment score, the element of the sliding window data structure is replaced by the element of the calculated alignment score. Therefore, after step 1, the sliding window data structure is "0, 5, 6, 9". Next, the sliding window slides along the calculation direction, and the sliding distance is the length of the sliding window data structure. At this point, the sliding window corresponds to "6, 0, 4". However, the last element does not correspond to an element within the calculation range. Therefore, the sliding window data structure should be compared with "6, 0, 4, #", where "#" represents an element outside the calculation range and is always less than the corresponding element in the sliding window data structure. Since the 0th element of the sliding window data structure is 0, and the corresponding element is 6, the 0th element should be updated. The other elements of the sliding window data structure remain unchanged. That is, after step 2, the sliding window data structure is "6, 5, 6, 9".

[0054] Then, a second maximum value operation (reduce_max, or dimensionality reduction and maximum value operation) is performed on the sliding window data structure to determine the maximum comparison score for this round. After step 2 above, the sliding window data structure is "6, 5, 6, 9". In step 3, the algorithm is used to determine the maximum value of these four elements in the sliding window data structure, which is "9".

[0055] Before round N, the maximum comparison score of the entire score matrix is ​​8. In round N, the maximum comparison score is 9, and the maximum comparison score of the entire score matrix needs to be updated to 9. In round N+1, the maximum comparison score is 5, and the maximum comparison score of the entire score matrix is ​​9, so it remains 9 without updating. In round N+2, the maximum comparison score is 6, and the maximum comparison score of the entire score matrix is ​​9, so it remains 9 without updating.

[0056] In some embodiments, the length of the sliding window data structure is equal to the data width of the SIMD instruction. Setting the length of the sliding window data structure to the data width of the SIMD instruction facilitates the execution of the SIMD instruction. Using SIMD instructions can improve computational efficiency. The two maximum value operations described above can use a single SIMD instruction (such as the SIMD MAX instruction) or a combination of SIMD instructions, which reduces the number of instructions used and improves computational efficiency compared to traditional sorting algorithms.

[0057] In some embodiments, after the calculation of the current round is completed, the comparison scores of the current round and the previous round are traversed along both ends of the anti-diagonal direction. If, at either end, the comparison scores of the current round and the previous round simultaneously show consecutive zero values, the calculation range is narrowed down to the non-zero score positions. The narrowed range is then compared with the calculation width parameter, and the smaller one is selected as the calculation range for the next round. This effectively prunes the calculation range for the next round based on the calculated comparison scores after each round of calculation, reducing the comparison scores required for the next round. Comparison scores continue to be calculated within the pruned calculation range, while in the pruned areas, no actual calculation is performed; instead, a value of 0 is directly assigned. While this method sacrifices some low-value details, it preserves high-value information (i.e., the optimal matching sequence), ensuring that the comparison results are not affected, while simultaneously reducing computational load and improving computational efficiency.

[0058] Specifically, after each round of calculation, the diagonal lines of the comparison scores calculated in this round and the diagonal lines of the comparison scores calculated in the previous round are traversed simultaneously from either end (i.e., from the bottom left corner or the top right corner). When both diagonal lines simultaneously show a continuous sequence of 0s starting from their respective starting points (i.e., both diagonal lines show several consecutive 0s starting from their respective starting points), it indicates that the portion of the comparison scores to be calculated in the next round belongs to a low-value region, and this region can be pruned to improve computational efficiency.

[0059] For example, the elements of the score matrix M corresponding to the comparison score calculated in the previous round are: The comparison scores calculated in this round correspond to the elements of the score matrix M as follows: ,and At this point, if we consider the elements from the starting point of the diagonal line representing the comparison scores calculated in the previous round to... All are 0, and the elements from the beginning of the slash of the comparison scores calculated in this round to... All values ​​are 0, and there are no non-zero values ​​separating the elements in the first region. Areas deemed to be of low value will be pruned in the next round of calculations.

[0060] The elements of the score matrix M corresponding to the comparison score calculated in the previous round are: The comparison scores calculated in this round correspond to the elements of the score matrix M as follows: ,and At this point, if we consider the elements from the starting point of the diagonal line representing the comparison scores calculated in the previous round to... All are 0, and the elements from the beginning of the slash of the comparison scores calculated in this round to... All values ​​are 0, and there are no non-zero values ​​separating the elements in the first region. Areas deemed to be of low value will be pruned in the next round of calculations.

[0061] Figure 5 illustrates a schematic diagram of the dynamic pruning process according to an embodiment of this application. Figure 5 provides an example of pruning two elements in the score matrix.

[0062] First, the traversal direction is from the bottom left corner to the top right corner. Because the first element of the two diagonals calculated in the third round ("0, 2, 2, 0") and the fourth round ("0, 4, 0")... , All are 0, therefore, correspondingly, the first element of the 5th round The branch is pruned. Then, the traversal direction is from the top right corner to the bottom left corner. In the third round of calculation, the diagonal lines "0, 2, 2, 0" and the fourth round of calculation, the first element of each diagonal line... , All of them are also 0, therefore, correspondingly, the first element of the 5th round The branches were pruned.

[0063] Therefore, the calculations in round 5 include: those that need to be actually calculated. and and directly assigning a value of 0 and After the fifth round of calculations is completed, based on the diagonal lines calculated in the fourth and fifth rounds, the elements for pruning in the sixth round are determined before the sixth round of calculations begins. This process is repeated for each round.

[0064] by For example, without pruning, if If an element is one of the elements in the optimal matching sequence of the query sequence obtained by backtracking, then after... Possible backtracking paths include , as well as Among them, due to Outside the scope of calculation, therefore It is not considered a possible backtracking path. The calculation yields 0, therefore and Both of these backtracking paths are This does not meet the requirements of a backtracking path. Therefore, a backtracking path must not include... Therefore, Pruning, which assigns a value of 0 directly without actual calculation, does not affect the final backtracked path (i.e., the optimal matching sequence).

[0065] In some embodiments, the above pruning method is used to narrow the calculation range of the next round. When the calculation range of the next round is narrowed to 0, the calculation of the next round and subsequent rounds is stopped. The comparison position corresponding to the maximum comparison score of the current record is used as the sum score matrix to determine the best matching sequence between the query sequence and the reference sequence.

[0066] Without pruning, multi-round computation ends when all rounds are completed. With pruning, it ends when all rounds are completed, or when the computation range for the next round is pruned (i.e., the actual computation range is reduced to 0). If the computation range for the next round is pruned, the diagonal regions in the next and subsequent rounds can be considered low-value regions, eliminating the need for the Smith-Waltman algorithm. Therefore, when later rounds are all low-value regions, pruning can avoid computation in these rounds, allowing for a faster completion of the score comparison process and improving computational efficiency in the seed expansion phase.

[0067] In some embodiments, for each element within the calculation range of the current round in the score matrix, the following formula is used for calculation:

[0068]

[0069] in, For the reference sequence The element and the query sequence of the first element The matching score of each element. Penalty score. Matching score. The score can be set according to actual needs. For example, when the two bases corresponding to the reference sequence and the query sequence are the same (e.g., AA), the matching score is a positive integer (e.g., +1), while when the two bases corresponding to the reference sequence and the query sequence are different (e.g., AC), the matching score is a negative integer (e.g., -1). In some embodiments, the matching score can be set as a substitution matrix, which is a 4x4 matrix. Each element in the matrix represents a matching score for a matching relationship, such as AT and TA, for which different matching scores can be set. Penalty score. This indicates the penalty for gaps. The penalty score can be further differentiated based on gap opening and gap extension, penalizing more discontinuous gaps versus fewer consecutive gaps. The penalty score can also be set according to actual needs. This application does not impose any restrictions on it.

[0070] In some embodiments, in step S150, the maximum alignment score of the score matrix is ​​traced back along the upper left, left, or upper direction to the element with an alignment score of 0. The starting point of this traversal path is the approximate maximum alignment score of the entire score matrix calculated above, and the traversal direction is the upper left, left, or upper direction (based on the approximate maximum alignment score of the entire score matrix calculated above). Starting from, respectively towards , and The direction of backtracking (the direction of the backtracking) terminates when the alignment score is 0. The backtracking direction can represent the relationship between two sequence bases: the top-left direction indicates a match (or replacement), the left direction indicates deletion, and the top direction indicates insertion.

[0071] In another embodiment of this application, an electronic device is also provided, which includes a memory and a processor. The memory stores computer-executable instructions that are executed by the processor. The computer-executable instructions include SIMD instructions. The processor executes the computer-executable instructions to implement the methods described in the embodiments of this application.

[0072] This electronic device employs an intra-query parallel computing strategy, which is suitable for scenarios with variable query sequences. It avoids waiting after the calculation of shorter query sequences. At the same time, it uses pruning strategies and heuristic search for extreme values ​​to reduce the actual calculation of elements in low-value regions and avoids extra traversal to determine the maximum value, thereby improving the computational efficiency of the seed expansion stage.

[0073] In another embodiment of this application, a computer-readable storage medium is also proposed, which stores a computer program or computer instructions that are executed by a processor to implement the technical solutions described in the embodiments of this application.

[0074] In another embodiment of this application, a computer program product is also proposed, which includes computer program code or computer instructions. When the computer program code or computer instructions are run, the technical solutions described in the embodiments of this application are implemented.

[0075] It will be understood by those skilled in the art that all or some of the steps and systems in the methods disclosed above can be implemented as software, firmware, hardware, and suitable combinations thereof. Some or all of the physical components can be implemented as software executed by a processor, such as a central processing unit, digital signal processor, or microprocessor, or as hardware, or as an integrated circuit, such as an application-specific integrated circuit. Such software can be distributed on a computer-readable storage medium. As is known to those skilled in the art, the term computer-readable storage medium includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information, such as computer-readable instructions, data structures, program modules, or other data. Computer-readable storage media include, but are not limited to, RAM, ROM, EEPROM, flash memory or other memory technologies, CD-ROM, digital versatile disc (DVD) or other optical disc storage, magnetic cartridges, magnetic tape, disk storage or other magnetic storage devices, or any other medium that can be used to store desired information and is accessible to a computer. Furthermore, it is known to those skilled in the art that communication media generally include computer-readable instructions, data structures, program modules, or other data in modulated data signals such as carrier waves or other transmission mechanisms, and may include any information delivery medium. Computer program products are software products that primarily implement their solutions through computer programs. The computer program code or computer instructions they contain can be obtained from computer-readable storage media and can also be transmitted, distributed, and downloaded over a network.

[0076] It should be noted that although the operations of the method of this application are described in a specific order in the accompanying drawings, this does not require or imply that these operations must be performed in that specific order, or that all the operations shown must be performed to achieve the desired result. On the contrary, the steps depicted in the flowchart can be performed in a different order. Additionally or alternatively, certain steps may be omitted, multiple steps may be combined into one step, and / or one step may be broken down into multiple steps.

[0077] It should be understood that when the terms "first," "second," "third," and "fourth," etc., are used in the claims, specification, and drawings of this application, they are used only to distinguish different objects and not to describe a specific order. The terms "comprising" and "including" as used in the specification and claims of this application indicate the presence of the described features, integrals, steps, operations, elements, and / or components, but do not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components, and / or collections thereof.

[0078] It should also be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the application. As used in this specification and claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used in this specification and claims refers to any combination and all possible combinations of one or more of the associated listed items, and includes such combinations.

[0079] Although the embodiments of this application are described above, the content is merely an example adopted for the purpose of facilitating understanding of this application and is not intended to limit the scope and application scenarios of this application. Any person skilled in the art described in this application may make any modifications and changes in the form and details of the implementation without departing from the spirit and scope disclosed in this application, but the scope of patent protection of this application shall still be determined by the scope defined in the appended claims.

Claims

1. A method for gene sequence alignment, comprising a seed search stage and a seed expansion stage, wherein, The seed expansion phase includes: obtaining multiple query sequences and their corresponding reference sequences to be compared based on multiple seeds determined for the read sequences to be compared via the seed lookup phase; expanding the length of these sequences to an integer multiple of the data width of the SIMD instruction and initializing the score matrix for each query sequence and its corresponding reference sequence; determining the calculation range of the current round in each score matrix, and calculating the comparison score of the query sequence and its corresponding reference sequence using the SIMD comparison instruction along the anti-diagonal direction of the score matrix; in response to the maximum comparison score of the current round being greater than the maximum comparison score of the previous round, obtaining the comparison position corresponding to the maximum comparison score of the current round; and after all rounds have been calculated, determining the best matching sequence between the query sequence and the reference sequence based on the comparison position corresponding to the maximum comparison score and the score matrix.

2. The method according to claim 1 further includes determining the maximum comparison score for the current round after the comparison score for the current round has been calculated, which includes: For the comparison scores already calculated in the current round, arranged along the anti-diagonal direction of the score matrix, a sliding window data structure is used to perform an element-by-element maximum value operation, and each maximum value of the operation is recorded in the sliding window data structure. The length and step size of the sliding window data structure are equal to the data width of the SIMD instruction. And move the sliding window data structure along the anti-diagonal direction to traverse all the comparison scores of the current round, thereby determining the maximum comparison score of the current round.

3. The method according to claim 2, wherein, The direction of moving the sliding window data structure is from the bottom left corner to the top right corner.

4. The method according to claim 1, wherein, Determining the calculation range for the current round includes determining the calculation range based on a pre-set calculation width parameter, which is set based on the matching position of the query sequence and the reference sequence.

5. The method according to claim 4, further comprising: After the calculation of the current round is completed, the comparison scores of the current round and the previous round are traversed along both ends of the anti-diagonal direction. If the comparison scores of the current round and the previous round both show consecutive zero values ​​at either end, the calculation range is narrowed down to the non-zero score position. The reduced range is compared with the calculated width parameter, and the smaller one is selected as the calculation range for the next round.

6. The method according to claim 5, further comprising: When the calculation range is reduced to 0, the calculation for the next round and subsequent rounds is stopped. The comparison position corresponding to the maximum comparison score of the current record is used as the sum score matrix to determine the best matching sequence between the query sequence and the reference sequence.

7. The method according to claim 1, wherein, After determining the calculation range for the current round in the score matrix, the SIMD read command is used to read the data containing the query sequence and reference sequence, which include the parts to be compared.

8. An electronic device, characterized in that, The method includes a memory and a processor, the memory storing computer-executable instructions that are executed by the processor, the computer-executable instructions including SIMD instructions, the processor executing the computer-executable instructions to perform the method according to any one of claims 1 to 7.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program or computer instructions that are executed by a processor to implement the method of any one of claims 1 to 7.

10. A computer program product, characterized in that, The computer program product includes computer program code or computer instructions, which, when executed, implement the method described in any one of claims 1 to 7.

Images