A method for improving gene sequence alignment based on BWT transformation by using multiplication algorithm

Abstract

The application relates to the field of bioinformatics, and relates to a method for improving gene sequence alignment based on BWT transformation by using a multiplication algorithm. First, a suffix array is directly calculated by using the multiplication algorithm through a reference long sequence, and then data structures required by BWT transformation are obtained. The BWT matrix does not need to be stored in space, N*N two-dimensional space is saved, and time is greatly improved. Seed sequences are subjected to pretreatment operation, the seed sequences are sorted, the same suffixes are gathered together, the intermediate result of the previous seed sequence alignment is conveniently saved, and the number of alignment times is reduced. After the sa interval is found, a search is conducted by using a BFS (Breadth-First Search) method, the BFS can guarantee that the optimal sequence score is found as soon as possible, the optimal sa is not missed, then a heap data structure is used for storage, the top of the heap is the highest score each time, the highest score item is obtained each time to conduct the next layer search, and the optimal matching result is obtained until the matching is finished.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to an improved gene sequence alignment method based on BWT transform using a doubling algorithm. Background Technology

[0002] Sequence alignment is a crucial and fundamental research area in bioinformatics, serving as the first step in gene assembly. It involves aligning multiple reads (seed sequences) obtained by a sequencer onto a reference genome. First-generation sequencing technologies, which acquired relatively small amounts of data, necessitated the design of software based on hash algorithms. However, with the rapid development of next-generation sequencing technologies, the ever-increasing volume of sequencing data presents a significant challenge to sequence alignment.

[0003] Therefore, the biological community has proposed using an algorithm based on the Burrows-Wheeler Transform (BWT) to solve sequence alignment, which can address the problem of large sequence space requirements and improve alignment speed while reducing space usage. However, the explosive growth of DNA sequencing data far exceeds its processing speed, and how to further improve alignment speed and reduce space usage remains a challenge in bioinformatics. Summary of the Invention

[0004] The purpose of this invention is to propose an improved gene sequence alignment method based on BWT transformation using a doubling algorithm, which further improves the sequence alignment speed and space reduction, and enhances the alignment accuracy.

[0005] To achieve the above objectives, the technical solution adopted in this invention is as follows: an improved gene sequence alignment method based on BWT transformation using a doubling algorithm, the main steps of which are:

[0006] Step 1: Calculate the sa suffix array for the reference long sequence using the doubling algorithm. The sa suffix array records the row number before sorting the N*N BWT matrix, i.e., the position of column F; where N is the length of the reference long sequence.

[0007] Step 2: Iterate through the long reference sequence Reference to create an array C. Calculate the number of elements smaller than character 'a' in the Reference[0,N-2] array C(a), and record the starting address of character 'a'.

[0008] Step 3: Calculate column L using the sa suffix array and the reference long sequence Reference, which is the rightmost column of the sorted N*N BWT matrix;

[0009] Step 4: Using the established C array, obtain column F, which is the leftmost column of the sorted N*N BWT matrix;

[0010] Step 5: Calculate the OCC array using column L from step 3. OCC(a,i) is the number of times 'a' appears before row i+1 in array L. Thus, the data structure required for BWT is completed.

[0011] Step 6: Sort the seed sequence of reads from right to left in lexicographical order using a combination of radix sort and counting sort, so that reads with the same suffix are grouped together;

[0012] Step 7: Find the same suffix by comparing the two reads above and below. Starting from the right, the first different base is the starting position of the alignment. After determining the starting position of the read alignment, use the BWT data structure constructed in Steps 1-5 and the BFS wide search method to find the sa interval, and output the matching position, the number of mismatches, and the running time.

[0013] This invention improves the gene sequence alignment method based on the Burrows-Wheeler Transform (BWT) by utilizing a doubling algorithm. It maps short sequences to long sequences, finding the optimal matching position of the short sequence on the long sequence. The doubling algorithm is used to improve the preprocessing stage of BWT. Unlike traditional BWT, which involves rotating to construct matrices and calculating various data structures, the improved method directly calculates the suffix array from the reference long sequence and quickly calculates the required data structure based on the suffix array. This not only eliminates the need to store the BWT matrix spatially, saving N*N (where N is the length of the reference gene sequence) two-dimensional space, but also significantly improves time efficiency, as it avoids the time-consuming operations of rotating and sorting the gene sequence as in the original BWT. Furthermore, this invention preprocesses the seed sequences, sorting them to group reads with the same suffix together. This facilitates saving intermediate results (with the same suffix) from the previous seed sequence alignment, reducing the number of alignments.

[0014] Compared with existing methods, the beneficial effects of the present invention are mainly reflected in:

[0015] (1) By using the doubling algorithm to build a suffix array sa for the long reference sequence, space and time complexity are saved because it is not necessary to build an N*N matrix (where N is the length of the reference) of the traditional BWT. If the traditional method is to build a BWT matrix and then sort the BWT matrix to obtain the sa suffix array, the time and space complexity will be too large because the length of the reference is too large. Sorting the N rows of the BWT matrix would have a time complexity of at least O(N*N*logN). However, the time complexity of directly sorting the suffix of the reference to obtain the sa suffix array using the doubling algorithm will be O(N*logN), and N*N space will be saved.

[0016] (2) After finding the sa interval, use the BFS breadth search method to search. Breadth search can ensure that the best sequence score is found as quickly as possible without missing the best sa. Then use a heap data structure to store it. Each time, the top of the heap will be the highest score. Each time the highest score item is obtained, the next level of search is performed until the matching ends, which is the best matching result.

[0017] (3) Group reads with the same suffix together and save the intermediate results of the comparison to reduce the number of accesses with low cache hit rate.

[0018] (4) Experiments have shown that the present invention improves the gene sequence alignment method based on BWT transformation by using a doubling algorithm. When the reference sequence is large, the doubling algorithm significantly improves the processing time for establishing the sa suffix array. At the same time, the present invention also has certain advantages over existing alignment methods and software in terms of both the time required for precise alignment and the accuracy of fuzzy alignment. Attached Figure Description

[0019] Figure 1 This is a flowchart of the algorithm of the present invention.

[0020] Figure 2 Using Reference: "ACAACG$" (where $ is the terminating character) as an example, create a BWT matrix, shift each row one character to the left, and then sort the BWT matrix.

[0021] Figure 3 For BWT-related data structures (using Reference: "ACAACG" as an example).

[0022] Figure 4 A comparison chart showing the reads before and after sorting.

[0023] Figure 5 This describes the matching process using "ACAACG" as an example and "AAC" as read. Detailed Implementation

[0024] Please see Figure 1 As shown, the gene sequence alignment method based on BWT transform proposed in this invention, which utilizes a doubling algorithm to improve the method, mainly comprises the following steps:

[0025] Step 1: Calculate the sa suffix array for the reference long sequence using the doubling algorithm. The sa suffix array records the row number before sorting the N*N BWT matrix, i.e., the position of column F; where N is the length of the reference long sequence.

[0026] Step 2: Iterate through the long reference sequence Reference to create an array C. Calculate the number of elements smaller than character 'a' in the Reference[0,N-2] array C(a), and record the starting address of character 'a'.

[0027] Step 3: Calculate column L using the sa suffix array and the reference long sequence Reference, which is the rightmost column of the sorted N*N BWT matrix;

[0028] Step 4: Using the established C array, obtain column F, which is the leftmost column of the sorted N*N BWT matrix;

[0029] Step 5: Calculate the OCC array using column L from Step 3. OCC(a,i) represents the number of times 'a' appears before row i+1 in array L. This completes the data structure required for BWT (BWT related data structures are as follows). Figure 3 (as shown);

[0030] Step 6: Sort the seed sequence of reads from right to left in lexicographical order using a combination of radix sort and counting sort, so that reads with the same suffix are grouped together;

[0031] Step 7: Find the same suffix by comparing the two reads above and below. Starting from the right, the first different base is the starting position of the alignment. After determining the starting position of the read alignment, use the BWT data structure constructed in Steps 1-5 and the BFS wide search method to find the sa interval, and output the matching position, the number of mismatches, and the running time.

[0032] The specific implementation algorithm of this invention is described below:

[0033] 1. Obtain the suffix array (sa) using the doubling algorithm.

[0034] Taking the sequence ACAACG as an example, its suffixes are ACAACG,CAACG,AACG,ACG,CG,G. The suffix array sa[i] is the index of the starting position of the suffix with rank i in the Reference. It actually records the actual position of the suffix with rank i in the Reference. Because the positions are disordered after sorting the BWT matrix, the suffix array is used to record the actual position.

[0035] The doubling algorithm used in this invention mainly utilizes two properties: First, the first k characters of suffixes i and j are equal, and the first k characters of suffixes i+k and j+k are equal. Second, the first 2k characters of suffix i are less than the first 2k characters of suffix j. Knowing the order of the first k characters of all suffixes, the order of the first 2k characters of all suffixes can be determined. The ranking of the first k characters of suffix i is used as the primary key, and the ranking of the first k characters of suffix i+k is used as the secondary key. Then, radix sort is used to calculate the ranking of the first 2k characters of suffix i.

[0036] The following uses the string Reference="ACAACG$" as an example to illustrate the detailed process of constructing a suffix array using a doubling algorithm ($ is smaller than any other character, omitted in the table, and replaced with 0 for sorting):

[0037] Table 1. Doubling Algorithm Process

[0038]

[0039] As can be seen from Table 1, in the first sorting, A has the smallest label of 1, G has the largest label of 3, and C has the next largest label of 2.

[0040] The first doubling increment is 1, which is the suffix formed by combining the characters in the first and second columns (AC = 12), the suffix formed by combining the characters in the second and third columns (CA = 21), and so on until the last column (G$ = 30). Then, radix sort is performed on these 6 two-digit numbers. It is not necessary to perform radix sort on all the numbers because the previous sorting already revealed the number of identical first keys (the tens digit of the two-digit number) and their sorting positions. Only the numbers with identical first keys need to be sorted. If there are x numbers with identical first keys, the subsequent sorted numbers are incremented by (x-1).

[0041] The second doubling increment is 2, which is the suffix of the first and third columns combined (acaa, i.e., 21), the suffix of the second and fourth columns combined (caac, i.e., 32), all the way to the last column (G$, i.e., 50), and then the radix is ​​sorted again (the third sort).

[0042] The sorting process ends when all the resulting numbers are distinct, because doubling the sorting process would result in the tens digit of the two-digit number being the same as the previous sorting result. For example, in Table 1, the third sorting result has no duplicate numbers, and doubling the sorting process again would result in the fourth sorting result being 2, 4, 1, 3, 5, 6.

[0043] Each time the sorting length is multiplied by 2, that is, it is continuously doubled, and it takes a maximum of log(N) times to complete the sorting, while the traditional sorting of the BWT matrix requires N (N is the length of the Reference) sorting times.

[0044] The resulting sorted array is denoted as the rank[i] array, where the index i represents the position of the key and the value represents the size of the key. Simply put, the sa suffix array indicates who is ranked which number, while rank indicates "your ranking".

[0045] Taking Reference="ACAACG$" as an example, after the doubling algorithm, we get rank[i] = 2, 4, 1, 3, 5, 6.

[0046] rank[0] = 2 indicates that the key of the second index of the sa suffix array is 0.

[0047] rank[1] = 4 indicates that the key of the 4th index of the sa suffix array is 1.

[0048] rank[2] = 1 indicates that the key of the first index of the sa suffix array is 2.

[0049] rank[3] = 3 indicates that the key of the third index of the sa suffix array is 3.

[0050] rank[4] = 5 indicates that the key of the 5th index of the sa suffix array is 4.

[0051] rank[5] = 6 indicates that the key of the 6th index of the sa suffix array is 5.

[0052] Depend on Figure 2 It can be seen that in the traditional matrix construction method, after left shifting, the last row of the matrix (row 6, starting from 0) must begin with the terminating character $. Since $ is the smallest, the first character of row 0 after sorting the matrix must be $, i.e., sa[0] = 6 (i.e., the length of Reference). From this, we can obtain sa[i] = 6, 2, 0, 3, 1, 4, 5. Therefore, it can be seen that the present invention can obtain the sa suffix array without using the traditional method of constructing a BWT matrix.

[0053] 2. Use the Reference and sa suffix arrays to obtain all the remaining data structures of the BWT (Burrows-Wheeler Transform).

[0054] The C(a) array calculates the number of characters smaller than a in the Reference sequence. The size is $<A<C<G<T$. The actual meaning of C(a) is to calculate the starting address of character a. C(a)+1 is the starting address of character a. Since C(A) calculates the number smaller than A, when encountering character A, the corresponding C(C), C(G), and C(T) for C, G, and T are all larger than C(A) because C, G, and T are all larger than A. So, the corresponding C(C), C(G), and C(T) for C, G, and T all need to be incremented by 1. The time complexity is O(N), where N is the length of the Reference.

[0055] The F array represents the first column after the BWT matrix is sorted. So, the same characters are grouped together. The F array can be directly obtained through the established C array above. F[0] = '$', F[1] = 'A', F[x] = 'C', where the value of x is C[c], F[y] = 'G', where the value of y is C[G]. Calculate the F array in sequence.

[0056] From Figure 3 the relationship in the BWT matrix, it can be seen that the F column and the L column must be adjacent in the original gene sequence. Since the sa suffix array represents the positions of the F array, and the L array is the previous character of the F array, only need to directly intercept the original gene sequence according to the positions of sa.

[0057] After the L array is established, the OCC array is established through the L array. OCC(a, i) is a two-dimensional array, which represents the number of times character a appears before the L column in the (i + 1)-th row. Taking Figure 3 OCC(A, 5) = 3 in it as an example, it means that it has appeared 3 times before the 6th row in the L column.

[0058] 3. Sort multiple reads seed sequences

[0059] The present invention uses an algorithm that combines radix sort and counting sort to sort the reads, and this algorithm has a high sorting efficiency for reads.

[0060] First step: Establish a hash index of the reads data. Read the reads file, generate forward and reverse reads, and count the number of occurrences of the same hash key value to establish a hash index of the reads data.

[0061] The second step involves loading the corresponding reads data blocks based on the reads hash index. Then, the reads data blocks are sorted using radix sort and counting sort algorithms. For reads within the same hash block (i.e., a group of reads with the same first Km), after excluding Km, they are divided into several segments of 8bp length, and radix sort is performed (e.g., if the read length is 48bp and K=16, the remaining 32bp is divided into 4 groups, each 8bp long). Counting sort is used for each 8bp segment.

[0062] During the counting sorting process, sometimes the data blocks are very large, such as millions of data points, and sometimes the data blocks are very small, with only a few dozen data points. This results in data that is sometimes dense and sometimes sparse. In the case of dense data, this invention adopts a method similar to Km hash classification sorting. When the data is sparse, the 8bp is divided into two groups of 4bp each, and radix sort is performed. When the data volume is less than 64bp, quicksort can be used directly.

[0063] The time complexity of sorting reads is O(L*K), where L is the length of the reads and K is the size of the reads data blocks. A characteristic of sorted reads is that reads with the same suffix will be grouped together, such as... Figure 4 As shown in the diagram. This clusters reads with high similarity together, reducing the number of comparisons required.

[0064] 4. Locate the sa interval and perform sequence alignment.

[0065] After the BWT-related data structures are established, the C array is used to find the start and end positions of the first round of the sa suffix array. `sp` represents the start position of the sa interval, and `ep` represents the end position of the sa interval. The first read sequence is compared from right to left starting from the last element, and then compared with the next read to obtain the same suffix intervals. The results of the same suffix alignments (sa intervals) are saved. For the alignment of the next read, its starting position is the length of the read minus the same suffix interval as the previous read; that is, the first different base from right to left is the starting position of the read alignment.

[0066] The first round of sp and ep (i represents the starting position of the seed sequence alignment):

[0067] sp = C[i]

[0068] ep = C[i+1]-1

[0069] Then, the sa interval is calculated using the backward search algorithm formula and stored using a heap data structure. The backward search algorithm formula is as follows:

[0070] sp(iQ)=C[i]+OCC[i,sp(Q)-1]

[0071] ep(iQ)=C[i]+OCC[i,ep(Q)]-1

[0072] Figure 5 To illustrate this matching process more clearly, let's take Reference="ACAACG" and read="AAC" as an example.

[0073] The backward search algorithm is used to find the matching interval *sa*. Within *sa*, a matching sequence is searched, and then a wide-first search (BFS) using a heap data structure is employed to find the difference where spaces can be inserted, allowing for mismatches. Utilizing the BFS method guarantees finding the optimal sequence score in the fastest possible time.

[0074] The steps of the breadth search are as follows:

[0075] Step 1: Initialize and create a priority queue substrHeap. The initial sequence is empty, meaning the initial score is 0.

[0076] Step 2: Initialize a Heap variable (initialize it to 0), store the results of each search, and score each result (the scoring standard is that the fewer the number of inserted spaces and mismatches, the higher the score), and the top of the heap is the item with the highest score.

[0077] Step 3: Calculate the `sa` interval using the backward search algorithm formula, then determine if this `read` is a substring of the `Reference`. If not, insert spaces or mismatches and try matching again. If the match exceeds the upper bound of the error, the match fails. If a match is successful within the error range, store the intermediate results in a heap and continue matching, using the highest-scoring item each time to proceed to the next level of search.

[0078] Step 4: Judge mismatch, gap open, match, etc., and score them. Then save the data in a heap data structure. The top of the heap will always be the highest score, which is the optimal matching score. The final output is also output using the data stored in the heap.

[0079] Example

[0080] The reference genome segment used in this example was selected from chromosome 1 of the Williams 82 soybean cultivar. The genome data file was downloaded from the official website of the National Center for Biotechnology Information (http: / / www.ncbi.nlm.nih.gov). The sequence number is >NC_016088.4Glycine max cultivar Williams 82chromosome 1,Glycine_max_v4.0,wholegenome shotgun sequence. The segment length is 50,000 bp. The comparison gene segment was selected from chromosome 1LG6 of the pea genome assembly. The sequence number is >LR535892.1Pisum sativum genome assembly,chromosome:1LG6,whole genome shotgun sequence, with each short sequence segment being 100 bp in length.

[0081] 1. Time comparison of whether the doubling algorithm is used to build the sa suffix array (only in the preprocessing stage of Reference, compared with traditional BWT):

[0082] Table 2 Comparison of whether doubling was used in time

[0083]

[0084] The results in Table 2 show that when the reference is large, the time for preprocessing (building the sa suffix array) by doubling algorithm is greatly improved, while the traditional BWT preprocessing time is very slow.

[0085] 2. Time comparison for precise matching:

[0086] Table 3. Time Comparison Results of Three Software Programs

[0087]

[0088] In terms of the comparison time of the three software methods, the present invention is significantly better than Bowtie2 and BWA software in terms of the comparison time of a lower number of reads. However, when the number of reads is increased, the results are not the best, and the results are lower than BWA software but higher than Bowtie2 software.

[0089] 3. Comparison of accuracy in fuzzy alignment:

[0090] Table 4. Accuracy Comparison Results of Three Software Programs

[0091]

[0092]

[0093] In terms of the comparison accuracy of the three methods, the present invention performs well overall, especially with high accuracy under low read conditions. However, as the number of reads increases, the comparison accuracy decreases slightly, indicating that although the BFS wide search method speeds up the search, it still loses some accuracy.

[0094] Overall, this invention improves the gene sequence alignment method based on BWT transform by utilizing a doubling algorithm. When the reference sequence is large, the doubling algorithm significantly improves the processing time for establishing the sa suffix array through preprocessing. Furthermore, this invention also offers advantages over existing alignment methods and software in terms of both precise alignment time and accuracy in fuzzy alignment.

[0095] The above description is merely an example and illustration of the concept of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described or use similar methods to replace them, as long as they do not deviate from the concept of the invention or exceed the scope defined in the claims, they should all fall within the protection scope of the present invention.

Claims

1. A gene sequence alignment method based on BWT transform improved by utilizing a doubling algorithm, characterized in that, The comparison method steps are as follows: Step 1: Use the doubling algorithm to calculate the sa suffix array for the reference long sequence Reference. The sa suffix array records the row numbers before sorting the traditional N*N BWT matrix, that is, the positions of column F; where N is the length of the reference long sequence Reference; The suffix array sa[i] is the index of the starting position of the suffix ranked i in Reference, actually recording the actual position of the suffix ranked i in Reference. Because the positions are chaotic after sorting the BWT matrix, the suffix array is used to record the actual positions; The first k characters of suffix i and suffix j are equal, and the first k characters of suffix i + k and suffix j + k are equal, and the first 2k characters of suffix i are less than the first 2k characters of suffix j. Given the size relationship of the first k characters of all suffixes, we know the size relationship of the first 2k characters of all suffixes. Take the rank of the first k characters of suffix i as the first keyword, take the rank of the first k characters of suffix i + k as the second keyword, and then use radix sort to calculate the rank of the first 2k characters of suffix i; Step 2: Traverse the reference long sequence Reference to establish the C array. The C(a) array calculates the number of elements smaller than character a in Reference[0, N - 2], recording the starting address of character a; The C(a) array calculates the number of characters smaller than a in the Reference sequence. The size is $<A<C<G<T$. The actual meaning of C(a) is to calculate the starting address of character a. C(a) + 1 is the starting address of character a; because C(A) calculates the number smaller than A, when encountering character A, the corresponding C(C), C(G), C(T) of C, G, T are all larger than C(A). Since C, G, T are all larger than A, the corresponding C(C), C(G), C(T) of C, G, T all need to be incremented by 1; the time complexity of this is also O(N), where N refers to the length of Reference; Step 3: Use the sa suffix array and the reference long sequence Reference to calculate column L, that is, the rightmost column after sorting the traditional N*N BWT matrix; Step 4: Obtain column F through the established C array, that is, the leftmost column after sorting the traditional N*N BWT matrix; The F array represents the first column after sorting the BWT matrix, so the same characters are grouped together. Through the established C array above, the F array can be directly obtained. F[0] = ‘$’, F[1] = ’A’, F[x] = ‘C’, where the value of x is C[c], F[y] =‘G’, where the value of y is C[G], and calculate the F array in turn; Column F and column L must be adjacent in the original gene sequence; since the sa suffix array represents the positions of the F array, and the L array is the previous character of the F array, just directly intercept the original gene sequence according to the positions of sa; Step 5: Calculate the OCC array using column L from step 3. OCC(a, i) is the number of times 'a' appears before row i+1 in array L. Thus, the data structure required for BWT is completed. Step 6: Sort the seed sequence of reads from right to left in lexicographical order using a combination of radix sort and counting sort, so that reads with the same suffix are grouped together; The first step is to build a hash index for the reads data: read the reads file, generate reads in both directions, count the number of occurrences of the same hash key value, and build a hash index for the reads data; The second step is to load the corresponding reads data blocks according to the reads hash index, and then sort the reads data blocks using radix sort and counting sort algorithms. For reads within the same hash block, i.e., a group of reads with the same first Km, after excluding Km, they are divided into several segments of 8bp in length, and radix sort is performed. For example, if the read length is 48bp and K=16, the remaining 32bp is divided into 4 groups, each with a length of 8bp. Counting sort is used each time 8bp is sorted. Step 7: Find the same suffix by comparing the two reads above and below. Starting from the right, the first different base is the starting position of the alignment. After determining the starting position of the read alignment, use the BWT data structure constructed in Steps 1-5 and the BFS wide search method to find the sa interval, and output the matching position, the number of mismatches and the running time. After the BWT-related data structure is established, the C array is used to find the start and end positions of the first round of the sa suffix array; sp represents the start position of the sa interval and ep represents the end position of the sa interval; the first read sequence is compared from right to left starting from the last position of the read, and compared with the next read to obtain the same suffix interval, and the result of the same suffix comparison is saved as sa interval; for the comparison of the next read, its comparison start position is the length of the read minus the same suffix interval with the previous read, that is, the first different base from right to left is the start position of the read comparison; In the first round, sp and ep, i represents the starting position of the seed sequence alignment: Then, the sa interval is calculated using the backward search algorithm formula and stored using a heap data structure; the backward search algorithm formula is as follows: The backward search algorithm is used to find the matching interval 'sa', and then the matching sequence is searched within the interval 'sa'. Finally, a widened BFS using a heap data structure is used to find the difference where spaces can be inserted and mismatches are allowed. The Breadth-First Search (BFS) method can guarantee finding the optimal sequence score in the fastest way. The steps of the BFS are as follows: Step 7-1: Initialize and create a priority queue substrHeap. The initial sequence is empty, meaning the initial score is 0. Step 7-2: Initialize a Heap variable, initially set to 0, to store the results of each search, and score each result. The scoring criteria is that the fewer the number of inserted spaces and mismatches, the higher the score, and the top of the heap is the item with the highest score. Step 7-3: Calculate the sa interval using the backward-search algorithm formula, and then determine whether this read is a substring of the Reference. If not, insert a space or mismatch and try matching again. If the error exceeds the upper bound, the matching fails. If the matching is successful within the error, store the intermediate results in a heap and continue matching. Each time, obtain the item with the highest score and proceed to the next level of search. Step 7-4: Judge mismatch, gap open, match, etc., and score them. Then save them in a heap data structure. The top of the heap will always be the highest score, which is the optimal matching score. The final output is also output using the data stored in the heap.

Images