Sequence alignment method and seed search method and device thereof

Abstract

The invention discloses a sequence alignment method and a seed search method and device thereof. The seed searching method comprises the steps that an index structure is constructed based on a reference genome sequence S, and the constructed index structure comprises an FM-Index index, a TBWT array, a TOCC array and a CNT array. And taking out the last base from the base fragments to be compared as an initial matching target, and locking the current matching region according to the CNT array and the initial matching target. And starting from the last but one basic group in the basic group fragment, matching two basic groups with S according to the current matching area and the constructed index structure in a reverse order every time until the number of the matched basic groups in the basic group fragment exceeds a preset threshold value, and taking the matched base sequence in the base fragment as a seed sequence. According to the method, the number of times of irregular memory access is reduced, so that the efficiency of a seed search task is improved, and a more efficient gene sequence comparison scheme is provided.

Description

Technical Field

[0001] This application relates to the fields of biological sequencing and gene data analysis, specifically to methods and equipment for sequence alignment. Background Technology

[0002] The statements in this section are merely to provide background information in relation to this application to aid in understanding it, and such background information does not necessarily constitute prior art.

[0003] With the accelerated industrialization of precision medicine, gene data analysis plays an increasingly important role. Currently, second-generation short-read data processing remains the dominant method for gene data analysis in clinical practice. Statistics show that over 70% of gene sequencing data originates from short-read sequences generated by the Illumina sequencing platform. Sequence alignment, as a fundamental step in the gene data analysis process that connects raw sequencing data with subsequent variant detection and gene annotation, is not only a core step determining the accuracy of the analysis results but also the most time-consuming. For each gene sequence fragment (read) generated by the sequencing platform, the current standard alignment process mainly consists of four consecutive steps: seed sequence matching, seed sequence tandem alignment, filtering low-quality sequence strands, and sequence strand extension and scoring to determine the final genomic localization result.

[0004] In the above process, the step of finding matching seed sequences has become one of the main bottlenecks restricting the performance improvement of the entire alignment process because it requires feature extraction and preliminary matching of massive short read sequences on a base-by-base basis. Summary of the Invention

[0005] The purpose of this application is to provide a new seed search method and a sequence alignment scheme based thereon, so as to improve the efficiency of seed search matching and alignment process.

[0006] The objective of this application is achieved through the following technical solution: According to a first aspect of this application, a seed search method for sequence alignment is provided, comprising: constructing an index structure based on a reference genome sequence S, wherein the constructed index structure includes an FM-Index index and a TBWT array, a TOCC array, and a CNT array; the TBWT array consists of base character pairs consisting of the last two bases of each row of the BWT matrix in the FM-Index index; the TOCC array is used to record the number of times a specified base character pair appears in the first k rows of the TBWT; and the CNT array is used to record the starting position of each type of base in the first column of the BWT matrix. The last base is taken from the base fragment to be aligned as an initial matching target, and the current matching region is locked according to the CNT array and the initial matching target. Starting from the second-to-last base in the base fragment, two bases are taken in reverse order each time and matched with S according to the current matching region and the constructed index structure until the number of matched bases in the base fragment exceeds a preset threshold, and the matched base sequence in the base fragment is used as a seed sequence.

[0007] Preferably, locking the current matching region based on the CNT array and the initial matching target includes: obtaining the starting position sp of the initial matching target in the first column of the BWT matrix based on the CNT array; obtaining the ending position ep of the initial matching target in the first column of the BWT matrix based on the CNT array; and locking the current matching region [sp, ep-1] based on the starting position sp and the ending position ep.

[0008] Preferably, starting from the second-to-last base in the base fragment, two bases are taken in reverse order each time to match with S based on the current matching region and the constructed index structure, including: using the second-to-last base in the base fragment as the first index b1; using the base preceding the first index b1 in the base fragment as the second index b2; updating the current matching region [sp, ep-1] based on the current matching region combined with the TOCC array, the CNT array, and the first index b1 and the second index b2; using the base preceding the second index as the new first index, and returning to execute the step of updating the second index.

[0009] Preferably, based on the current matching region combined with the TOCC array, CNT array, and the first index b1 and the second index b2, the current matching region [sp, ep-1] is updated according to the following formula: sp =CNT[b2]+TOCC(sp-1,b1,b2), ep =CNT[b2]+TOCC(ep-1,b1,b2); Where CNT[b2] represents the starting position of b2 in the first column of the BWT matrix; TOCC(sp-1,b1,b2) represents the number of times the base pair formed by b1 and b2 appears in the first sp-1 rows of the BWT matrix; TOCC(ep-1,b1,b2) represents the number of times the base pair formed by b1 and b2 appears in the first ep-1 rows of the BWT matrix; and the region [sp ,ep -1] is used as the new current matching region [sp,ep-1].

[0010] Preferably, it further includes: mapping the current matching region to the corresponding position in the reference genome sequence based on the sampling suffix array of the FM-Index index, thereby obtaining the seed sequence and its starting position in the reference genome sequence.

[0011] Preferably, constructing an index structure based on a reference genome sequence S includes: performing a BWT transformation on the reference genome sequence S to obtain a BWT matrix corresponding to S; forming a BWT array from the last base of each row of the BWT matrix; forming a TBWT array from the base character pairs formed by the last two bases of each row of the BWT matrix; setting an OCC array to record the number of times a specified base appears in the first k rows of the BWT array; setting a TOCC array to record the number of times a specified base character pair appears in the first k rows of the TBWT array; setting a CNT array to record the starting position of each type of base in the first column of the BWT matrix; setting an SA array to record the corresponding position of each row of the BWT matrix in S; and forming an index structure corresponding to S from the BWT array, OCC array, Count array, SA array, TBWT array, and TOCC array.

[0012] Preferably, it further includes: taking the base sequence composed of the last N bases of each row of the BWT matrix to form an N-order BWT array, where N is a positive integer greater than 2; setting an N-order OCC array, which is used to record the occurrence count of a specified base sequence in the first k rows of the N-order BWT array.

[0013] Preferably, the method further includes compressing and storing the index structure in the following manner: setting one or more checkpoints and the interval M between checkpoints; traversing the rows of the BWT matrix, if the row number of the current row is an integer multiple of the interval M, triggering checkpoint storage, and storing the data of the OCC array, TOCC array, and TBWT array corresponding to the current row into a checkpoint; setting a MOCC array in the checkpoint to store the occurrence count of all base character pairs recorded in the TOCC array within a preset range, wherein the preset range includes each row from the row where the checkpoint is located to the row between the checkpoint and the next checkpoint; and storing all checkpoints.

[0014] Preferably, the method further includes obtaining TOCC data from the compressed index structure in the following manner: calculating the position of the checkpoint corresponding to the current row based on the row number of the current row and the interval M between checkpoints; retrieving the value corresponding to the specified base character pair from the TOCC array stored at the corresponding checkpoint based on the position of the checkpoint, and assigning it to tocc; retrieving the value corresponding to the specified base character pair from the MOCC array corresponding to the current row, summing tocc with the value to obtain a new tocc value; and using the new tocc value as the TOCC data of the current row.

[0015] According to a second aspect of this application, a sequence alignment method is provided, comprising: determining the position of a seed sequence of a base fragment to be aligned in a reference genome using the seed search method in the first aspect; based on the position of the seed sequence in the reference genome, extending the seed sequence to both sides or one side and performing base-by-base alignment with the corresponding region of the reference genome sequence S, and calculating the alignment score of the extended region of each seed sequence based on the alignment results; evaluating the scores of all extended regions, and taking the match with the highest score as the output result.

[0016] According to a third aspect of this application, a computer device is provided, including a memory and a processor, wherein a computer program capable of running on the processor is stored in the memory, wherein the processor executes the program to implement the methods as described in the first and second aspects.

[0017] According to a fourth aspect of this application, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores a computer program that, when executed by a processor, implements the methods as described in the first and second aspects.

[0018] Compared with existing technologies, the main advantages of this application are: it enables a method to simultaneously expand two bases in a single memory read, thereby reducing the number of irregular memory accesses by half. By improving the shortcomings of irregular memory access, the efficiency of the seed search task in the gene sequence alignment process is improved, providing a more efficient gene sequence alignment scheme. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1 This is a schematic diagram of the data structure of the FM-Index index; Figure 2This is a schematic diagram of the data structure of the FMT-Index index according to an embodiment of this application; Figure 3 This is a schematic diagram of a checkpoint structure for compressed storage of an FMT-Index index according to an embodiment of this application; Figure 4 This is a flowchart illustrating a seed search method for sequence alignment based on the FMT-Index index according to an embodiment of this application. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided through specific embodiments. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of this application.

[0021] As mentioned in the background section, most mainstream short sequence alignment software (such as bwa-mem and bowtie2) uses algorithms based on the FM-Index index structure to achieve precise seed sequence matching. The FM-Index algorithm determines that when the number of consecutive matching bases in a short read sequence exceeds a preset threshold (usually set to 19), the continuous sequence is considered a valid seed sequence, which is then used for subsequent sequence chain extension and full sequence alignment.

[0022] The inventors discovered that the FM-Index algorithm is inherently a memory-intensive computational model, involving numerous irregular memory access operations during the matching process. These irregular memory accesses are characterized by: 1) random and unpredictable distribution of accessed memory addresses (for example, in the FM-Index algorithm, addresses calculated from the Count and OCC arrays may be scattered across various memory locations); 2) unpredictable address distribution, rendering the prefetcher ineffective; and 3) the ability to match only one base at a time, resulting in a large number of accesses despite only a small amount of data read / written per transaction. Each transmission requires address establishment and response waiting, severely diluting the effective data transfer rate of the bus. Numerous scattered memory access requests consume bus bandwidth, forcing subsequent requests to queue and lengthening the overall algorithm runtime. This directly leads to low algorithm efficiency, making it difficult to meet the high-throughput demands of gene data analysis. Therefore, addressing the drawbacks of irregular memory access is crucial for improving the efficiency of the seeding stage.

[0023] After in-depth research, the inventors believe that the key to the above problems lies in the inherent defects of the FM-Index data structure. Figure 1 This is a schematic diagram of the data structure for the FM-Index index. (Example:) Figure 1As shown, the FM-Index primarily comprises the BWT array, OCC array, Count array, and SA array. The BWT array consists of the last base of each row of the BWT matrix obtained by performing a BWT transformation on the reference genome sequence S; the OCC array records the frequency of a specified base in the first k rows of the BWT array; the CNT array records the starting position of each base class in the first column of the BWT matrix; and the SA array records the corresponding position of each row of the BWT matrix within S. It can be seen that the FM-Index records the relevant information of a single extracted base in the BWT matrix. Therefore, when using the FM-Index for seed sequence search and matching, only one base can be matched per memory read, resulting in very low matching efficiency.

[0024] To reduce the number of irregular accesses, this application proposes an improved index, FMT-Index. Figure 2 This is a schematic diagram of the data structure of the FMT-Index index according to an embodiment of this application. Figure 2 As shown, the FMT-Index data structure adds a TBWT array and a TOCC array to the FM-Index data structure. The TBWT array consists of base character pairs consisting of the last two bases of each row of the BWT matrix; the TOCC array records the number of times a specified base character pair appears in the first k rows of the TBWT array. In some embodiments, the TOCC array uses base character pair encoding to represent the number of times a specified base character pair appears in the first k rows of the TBWT array. For example, a single base can be encoded as two bits (A=00, C=01, G=10, T=11), so two consecutive bases can be encoded as 4 bits. The encoding of the first base is shifted left by 2 bits and then concatenated with the encoding of the next base.

[0025] According to one embodiment of this application, a method for constructing an FMT-Index index based on a reference genome sequence S is proposed. The method includes: First, performing a BWT transformation on the reference genome sequence S to obtain a BWT matrix corresponding to S. Then, forming a BWT array from the last base of each row of the BWT matrix, and forming a TBWT array from the last two base pairs of each row of the BWT matrix. An OCC array is set to record the frequency of a specified base in the first k rows of the BWT array. A TOCC array is set to record the frequency of a specified base pair in the first k rows of the TBWT array. A CNT array is set to record the starting position of each base type in the first column of the BWT matrix. An SA array is set to record the corresponding position of each row of the BWT matrix in S. Finally, the BWT array, OCC array, Count array, SA array, TBWT array, and TOCC array are combined to form the data structure of the FMT-Index index corresponding to S.

[0026] If the TOCC data is stored directly, the FMT-Index data structure will occupy a very large amount of memory (estimated at 480GB). Figure 3 This is a schematic diagram of a checkpoint structure for compressed storage of an FMT-Index index according to an embodiment of this application. (Refer to...) Figure 3 In some embodiments, a checkpoint compression strategy is used to store data in the FMT-Index. Specifically, this includes: setting one or more checkpoints and an interval M between them, where M is a positive integer; traversing the rows of the BWT matrix, and if the row number of the current row is an integer multiple of the interval M, triggering checkpoint storage, and storing the data from the OCC array, TOCC array, and TBWT array corresponding to the current row into a checkpoint; setting a MOCC array in the checkpoint to store the occurrence counts of all base character pairs recorded in the TOCC array within a preset range, where the preset range includes the row where the checkpoint is located to every row between the checkpoint and the next checkpoint; in some embodiments, the MOCC array is also used to store the occurrence counts of all bases recorded in the OCC array within a preset range; finally, all checkpoints are stored.

[0027] In some embodiments, the method further includes retrieving TOCC data from the compressed FMT-Index index as follows: Calculate the position of the checkpoint corresponding to the current row based on the row number of the current row and the interval M between checkpoints. Retrieve the value corresponding to the specified base character pair from the TOCC array stored for the corresponding checkpoint based on the checkpoint position, and assign it to tocc. Retrieve the value corresponding to the specified base character pair from the MOCC array corresponding to the current row, and sum tocc with this value to obtain a new tocc value. Use the new tocc value as the TOCC data for the current row.

[0028] Figure 4 This is a flowchart illustrating a seed search method for sequence alignment based on the FMT-Index index, according to an embodiment of this application. Figure 4 As shown, the method includes: Step S401, constructing an FMT-Index index based on the reference genome sequence S. Step S402, extracting the last base from the base fragment to be aligned as the initial matching target, and locking the current matching region based on the CNT array and the initial matching target. Step S403, starting from the second-to-last base in the base fragment, taking two bases in reverse order each time, and matching them with S based on the current matching region, the TOCC array, and the CNT array, until the number of matched bases in the base fragment exceeds a preset threshold, and using the matched base sequence in the base fragment as the seed sequence. This embodiment of the application implements a method that allows for simultaneous expansion of two bases in a single memory read based on the FMT-Index data structure, thereby reducing the number of irregular memory accesses by half. By improving the shortcomings of irregular memory access, the efficiency of the seed search task in the gene sequence alignment process is improved, providing a more efficient gene sequence alignment scheme.

[0029] In some embodiments, step S402, locking the current matching region based on the CNT array and the initial matching target b, may include: obtaining the starting position sp=CNT[b] of the initial matching target b in the first column of the BWT matrix based on the CNT array, and obtaining the ending position ep=CNT[b+1] of the initial matching target b in the first column of the BWT matrix based on the CNT array. Here, the base types include four types: A, C, G, and T. Base b+1 represents the base type that follows the base type of base b in the ACGT sorting. For example, if base b is A, then base b+1 is C; if base b is C, then base b+1 is G. Finally, the current matching region [sp, ep-1] is locked based on the starting position sp and the ending position ep.

[0030] In some embodiments, step S403 may include: using the penultimate base in the base fragment as a first index b1; using the base preceding the first index b1 in the base fragment as a second index b2; updating the current matching region [sp, ep-1] based on the current matching region combined with the TOCC array, the CNT array, and the first index b1 and the second index b2; using the base preceding the second index as a new first index, and returning to the step of updating the second index, until the number of matched bases in the base fragment exceeds a preset threshold (typically 19), and using the matched base sequence in the base fragment as a seed sequence.

[0031] In some embodiments, the current matching region [sp, ep-1] is updated based on the current matching region in combination with the TOCC array, the CNT array, and the first index b1 and the second index b2 according to the following formula: sp =CNT[b2]+TOCC(sp-1,b1,b2), ep =CNT[b2]+TOCC(ep-1,b1,b2); Where CNT[b2] represents the starting position of b2 in the first column of the BWT matrix; TOCC(sp-1,b1,b2) represents the number of times the base pair formed by b1 and b2 appears in the first sp-1 rows of the BWT matrix; TOCC(ep-1,b1,b2) represents the number of times the base pair formed by b1 and b2 appears in the first ep-1 rows of the BWT matrix; and the region [sp ,ep -1] is used as the new current matching region [sp,ep-1].

[0032] Because the FMT-Index data structure still includes an OCC array, and it resides in contiguous memory with the TOCC array, in some embodiments, both OCC and TOCC data can be retrieved simultaneously in a single memory access to match one or two bases. Specifically, the CNT and OCC arrays can be used to match one base at a time and update the matching region according to the following formula: sp =CNT[b]+OCC(sp-1,b), ep =CNT[b]+OCC(ep-1,b); Where CNT[b] represents the starting position of b in the first column of the BWT matrix; OCC(sp-1,b) represents the number of times the base b appears in the first sp-1 rows of the BWT matrix; OCC(ep-1,b) represents the number of times the base b appears in the first ep-1 rows of the BWT matrix.

[0033] In some embodiments, the method may further include: retrieving the corresponding value from the sampled suffix SA array based on the current matching region [sp, ep-1], forming a position array pos_arr, and using the position array pos_arr as the starting position of the seed sequence in the reference genome sequence S. This embodiment of the application finds the seed sequence and determines its starting position in the reference genome sequence S, thus completing the seeding stage.

[0034] To further reduce the number of irregular accesses, in some embodiments, the FMT-Index data structure may also include an N-order BWT array, composed of the last N bases of each row of the BWT matrix, where N is a positive integer greater than 2. Simultaneously, an N-order OCC array is set to record the occurrence count of a specified base sequence in the first k rows of the N-order BWT array. For example, if N=4, the FMT-Index data structure includes TBWT, 3BWT, and 4BWT, along with their corresponding TOCC, 3OCC, and 4OCC, all stored in contiguous memory. When performing seed search using the N-order OCC array, this data can be retrieved with a single memory access, matching 1 to N bases. The number of bases matched at one time can be flexibly selected, thereby further improving the efficiency and performance of the seed search task in the gene sequence alignment process.

[0035] According to another embodiment of this application, a sequence alignment method is provided, comprising: determining the position of a seed sequence of a base fragment to be aligned in a reference genome using the seed search method described in the foregoing embodiments; based on the position of the seed sequence in the reference genome, extending the seed sequence laterally or unilaterally and performing base-by-base alignment with the corresponding region of the reference genome sequence S; and calculating the alignment score of the extended region of each seed sequence based on the alignment results; evaluating the scores of all extended regions, and taking the match with the highest score as the output result.

[0036] According to yet another embodiment of this application, a computer device is provided, including a memory and a processor, wherein a computer program capable of running on the processor is stored in the memory, and wherein the processor executes the program to implement the method as described in the foregoing embodiments.

[0037] References to "various embodiments," "some embodiments," "one embodiment," or "embodiment," etc., in this specification refer to a specific feature, structure, or property described in connection with an embodiment that is included in at least one embodiment. Therefore, the appearance of the phrases "in various embodiments," "in some embodiments," "in one embodiment," or "in an embodiment," etc., throughout this specification does not necessarily refer to the same embodiment. Furthermore, specific features, structures, or properties can be combined in any suitable manner in one or more embodiments. Therefore, a specific feature, structure, or property shown or described in connection with one embodiment can be combined, in whole or in part, with features, structures, or properties of one or more other embodiments without limitation, provided that such combination is not illogical or inoperable.

[0038] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments or implementation methods. Similar or identical parts between the various embodiments of the present invention can be referred to mutually. The implementation principles and technical effects of the inventive concept can be referred to mutually, and will not be repeated here. Where there is no conflict, the various embodiments or implementation methods in this invention can be combined with each other.

[0039] It should be noted that although the steps are described in a specific order above, it does not mean that the steps must be executed in the above specific order. In fact, some of these steps can be executed concurrently, or even in a different order, as long as the required function can be achieved.

[0040] Embodiments of this application may be systems, methods, and / or computer program products. A computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for causing a processor to implement various aspects of this application.

[0041] Computer-readable storage media can be tangible devices that hold and store instructions for use by an instruction execution device. Computer-readable storage media can include, for example, but not limited to, electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination thereof.

[0042] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used in the embodiments of the present invention is chosen to best explain the principles, practical applications, or technological improvements in the market, or to enable other those skilled in the art to understand the disclosed embodiments of the present invention.

Claims

1. A seed search method for sequence alignment, wherein, include: An index structure is constructed based on the reference genome sequence S. The constructed index structure includes an FM-Index index and TBWT, TOCC, and CNT arrays. The TBWT array consists of base character pairs consisting of the last two bases of each row of the BWT matrix in the FM-Index index. The TOCC array is used to record the number of times a specified base character pair appears in the first k rows of the TBWT array. The CNT array is used to record the starting position of each type of base in the first column of the BWT matrix. The last base is taken from the base fragment to be compared as the initial matching target, and the current matching region is locked according to the CNT array and the initial matching target; Starting from the second-to-last base in the base fragment, two bases are taken in reverse order each time and matched with S according to the current matching region and the constructed index structure until the number of bases that have completed matching in the base fragment exceeds a preset threshold. The base sequence that has completed matching in the base fragment is then used as the seed sequence.

2. The method according to claim 1, wherein, Based on the CNT array and the initial matching target, the current matching region is locked, including: The starting position sp of the initial matching target in the first column of the BWT matrix is ​​obtained based on the CNT array; The end position ep of the initial matching target in the first column of the BWT matrix is ​​obtained based on the CNT array; The current matching region [sp, ep-1] is locked based on the starting position sp and the ending position ep.

3. The method according to claim 1, wherein, Starting from the second-to-last base in the base fragment, two bases are taken in reverse order each time and matched with S according to the current matching region and the constructed index structure, including: The penultimate base in the base fragment is used as the first index b1; The base preceding the first index b1 in the base fragment is used as the second index b2; Update the current matching region [sp, ep-1] based on the current matching region, combined with the TOCC array, CNT array, and the first index b1 and the second index b2; Use the base preceding the second index as the new first index, and return to perform the step of updating the second index.

4. The method according to claim 3, wherein, Based on the current matching region, combined with the TOCC array, CNT array, and first index b1 and second index b2, the current matching region [sp, ep-1] is updated according to the following formula: sp =CNT[b2]+TOCC(sp-1,b1,b2), Hey =CNT[b2]+TOCC(ep-1,b1,b2); Where CNT[b2] represents the starting position of b2 in the first column of the BWT matrix; TOCC(sp-1,b1,b2) represents the number of times the base pair consisting of b1 and b2 appears in the first sp-1 rows of the BWT matrix. TOCC(ep-1,b1,b2) represents the number of times the base pair consisting of b1 and b2 appears in the first ep-1 rows of the BWT matrix. The region [sp ,ep -1] is used as the new current matching region [sp,ep-1].

5. The method according to claim 1, further comprising: The current matching region is mapped to the corresponding position in the reference genome sequence based on the sampled suffix array of the FM-Index, thereby obtaining the seed sequence and its starting position in the reference genome sequence.

6. The method according to claim 1, wherein, The index structure constructed based on the reference genome sequence S includes: The reference genome sequence S is subjected to BWT transformation to obtain the BWT matrix corresponding to S; The last base of each row of the BWT matrix is ​​used to form a BWT array; The base pairs formed by the last two bases of each row of the BWT matrix are used to form the TBWT array; Set up the OCC array, which is used to record the number of times a specified base appears in the first k rows of the BWT array; Set up the TOCC array, which is used to record the number of times a specified base character pair appears in the first k rows of the TBWT array; Set up a CNT array, which is used to record the starting position of each type of base in the first column of the BWT matrix; Set up an array SA, which is used to record the corresponding position of each row of the BWT matrix in S; Combine the BWT array, OCC array, Count array, SA array, TBWT array, and TOCC array to form the index structure corresponding to S.

7. The method of claim 6, further comprising: Take the last N bases of each row of the BWT matrix to form a base sequence, and form an N-order BWT array, where N is a positive integer greater than 2; Set up an N-order OCC array, which is used to record the number of times a specified base sequence appears in the first k rows of the N-order BWT array.

8. The method of claim 6 further comprises compressing and storing the index structure in the following manner: Set one or more checkpoints and the interval M between checkpoints; Traverse the rows of the BWT matrix. If the row number of the current row is an integer multiple of the interval M, trigger checkpoint storage and store the data of the OCC array, TOCC array and TBWT array corresponding to the current row into a checkpoint. A MOCC array is set in the checkpoint to store the occurrence count of all base character pairs recorded in the TOCC array within a preset range, wherein the preset range includes the row where the checkpoint is located to each row between the checkpoint and the next checkpoint; Store all checkpoints.

9. The method of claim 8, further comprising retrieving TOCC data from the compressed index structure in the following manner: Calculate the position of the checkpoint corresponding to the current row based on the row number of the current row and the interval M between checkpoints; Based on the position of the checkpoint, retrieve the value corresponding to the specified base character pair from the TOCC array stored at the corresponding checkpoint and assign it to tocc; Extract the value corresponding to the specified base character pair from the MOCC array corresponding to the current line, and sum the tocc value with the value to obtain a new tocc value; Use the new TOCC value as the TOCC data for the current row.

10. A sequence alignment method, comprising: The seed sequence of the base fragment to be compared is determined in the reference genome using the seed search method described in claims 1-9; Based on the position of the seed sequence in the reference genome, the seed sequence is extended to both sides or one side and compared with the corresponding region of the reference genome sequence S base by base, and the alignment score of the extended region of each seed sequence is calculated based on the alignment results. Evaluate the scores of all extended regions and output the match with the highest score.

11. A computer device comprising a memory and a processor, wherein a computer program capable of running on the processor is stored in the memory, wherein the processor, when executing the program, implements the method according to any one of claims 1-10.

Images