Syncmer-based genome double-sequence alignment method and system

Abstract

The invention belongs to the technical field of gene sequence alignment, and particularly relates to a syncmer-based genome double-sequence alignment method and system.The method comprises the steps that sampling and subsequence extraction are conducted on the basis of a syncmer strategy, division is conducted according to the subsequences, a continuous long alignment interval is divided into a plurality of short alignment intervals, and the calculation efficiency is remarkably improved while high alignment accuracy is guaranteed; short comparison intervals are dynamically optimized and merged, so that the interval continuity and the comparison efficiency are effectively balanced; in order to reduce the influence of interval division on the comparison result, the local optimal comparison paths of the adjacent optimization comparison intervals are spliced according to the positions of the matched anchor points, so that a more accurate comparison result is generated.

Description

Technical Field

[0001] The present invention belongs to the technical field of gene sequence alignment, and specifically relates to a genomic double-sequence alignment method and system based on syncmer. Background Art

[0002] Sequence alignment is a fundamental task in bioinformatics, aiming to reveal evolutionary relationships, locate functional elements, and identify genetic variations by comparing the similarities of different sequences. Among them, double-sequence alignment, as the most fundamental and computationally intensive part of sequence alignment, directly affects the accuracy and efficiency of downstream tasks such as variant detection and gene annotation. How to break through its computational bottleneck and achieve efficient processing of large-scale genomic data is the core challenge in current bioinformatics research. With the progress of sequencing technology, the sequence length has leaped from kilobases to megabases, and traditional methods are difficult to apply due to excessive memory and long time consumption, and there is an urgent need for efficient algorithm innovation. Summary of the Invention

[0003] The present invention provides a genomic double-sequence alignment method and system based on syncmer.

[0004] The technical solution of the present invention is as follows: The present invention provides a genomic double-sequence alignment method based on syncmer, including: S1: Obtain a target sequence and a query sequence, and use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively. If the distance between two adjacent subsequences is greater than the distance threshold, insert at least one subsequence into the interval between the two adjacent subsequences so that the distance between the two adjacent subsequences is not greater than the distance threshold, to obtain the subsequence of the target sequence and the subsequence of the query sequence; After constructing hash indexes based on the subsequence of the target sequence and the subsequence of the query sequence respectively, obtain the target sequence index and the query sequence index; S2: According to the target sequence index and the query sequence index, match each subsequence of the query sequence in the target sequence, and record the two matching subsequences as matching anchors; Based on the matching anchors, divide the intervals of the query sequence and the target sequence respectively to obtain the alignment interval of the target sequence and the alignment interval of the query sequence; Optimize the alignment interval of the target sequence and the alignment interval of the query sequence respectively according to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlap situation between the alignment intervals, to obtain the corresponding optimized alignment intervals; S3: Each optimized alignment interval is processed by a double-sequence alignment algorithm to obtain a locally optimal alignment path; According to the positions of the matching anchors, splice the locally optimal alignment paths of adjacent optimized alignment intervals to generate an alignment result.

[0005] S2, according to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlap situation between alignment intervals, respectively optimize the alignment intervals of the target sequence and the alignment intervals of the query sequence to obtain corresponding optimized alignment intervals, specifically: For the target sequence or the query sequence, if the distance between two adjacent alignment intervals is less than a preset threshold, then merge the two adjacent alignment intervals to obtain the corresponding optimized alignment interval; wherein, the preset threshold is set according to the sequence characteristics of the alignment interval to be optimized; If the length of the alignment interval is less than the preset length, then merge this alignment interval with the adjacent alignment interval to obtain the corresponding optimized alignment interval; If there is an overlap between alignment intervals, then merge the overlapping alignment intervals to obtain the corresponding optimized alignment interval.

[0006] S3, according to the matching anchor position, splice the locally optimal alignment paths of adjacent optimized alignment intervals to generate an alignment result, specifically: According to the matching anchor position, if the offsets of the locally optimal alignment paths of adjacent optimized alignment intervals on the target sequence and the query sequence are the same, then merge the locally optimal alignment paths of adjacent optimized alignment intervals into a continuous alignment path; If there is an overlap or contradiction in the locally optimal alignment paths of adjacent optimized alignment intervals, then adopt a dynamic programming backtracking strategy to generate a globally optimal path based on the weights of the locally optimal alignment paths; After splicing the locally optimal alignment paths of adjacent optimized alignment intervals corresponding to all matching anchor positions, generate an alignment result.

[0007] S2, based on the matching anchors, respectively divide the query sequence and the target sequence into intervals to obtain the alignment intervals of the target sequence and the alignment intervals of the query sequence, specifically: Sort the matching anchors of the query sequence and the target sequence respectively. Taking two adjacent anchors as boundaries, after dividing the query sequence or the target sequence, take the termination position of the former of the two adjacent anchors as the starting position of the query sequence or the target sequence, and take the starting position of the latter of the two adjacent anchors as the termination position of the query sequence or the target sequence to obtain the alignment intervals of the target sequence and the alignment intervals of the query sequence.

[0008] S2, according to the target sequence index and the query sequence index, match each subsequence of the query sequence in the target sequence, and record the two successfully matched subsequences as matching anchors, specifically: Using the base sequence of the subsequence of the query sequence as the key, search in the hash index of the target sequence. If a matching subsequence is found, record the two successfully matched subsequences as matching anchors, and record the base sequence of the subsequence, the position of the subsequence in the target sequence, and the position of the subsequence in the query sequence.

[0009] In S1, syncmers are used to extract subsequences of the target sequence and the query sequence respectively, specifically: Select a k-mer within a sliding window containing multiple k-mers. If the s-mer is at the starting position or the end position of the k-mer, consider the k-mer as a subsequence; continue to slide the current window to extract subsequences until the subsequences of the target sequence and the query sequence are completely extracted.

[0010] In S1, based on the subsequences of the target sequence and the query sequence, hash indexes are constructed respectively. The base sequence of the subsequence of the target sequence is used as the key of the hash table, and the position where the subsequence of the target sequence appears in the target sequence is used as the value of the hash table; The base sequence of the subsequence of the query sequence is used as the key of the hash table, and the position where the subsequence of the query sequence appears in the query sequence is used as the value of the hash table.

[0011] The alignment result includes the global alignment score, the global alignment path, the alignment interval, the number of matching bases, the number of mismatched bases, and the number of inserted or deleted bases.

[0012] The present invention also provides a syncmer-based genomic dual-sequence alignment system, including: Index construction module: used to obtain the target sequence and the query sequence, and use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively. If the distance between two adjacent subsequences is greater than the distance threshold, insert at least one subsequence in the interval between the two adjacent subsequences so that the distance between the two adjacent subsequences is not greater than the distance threshold, and obtain the subsequences of the target sequence and the query sequence; After constructing hash indexes based on the subsequences of the target sequence and the query sequence respectively, obtain the target sequence index and the query sequence index; Alignment interval optimization module: according to the target sequence index and the query sequence index, match each subsequence of the query sequence in the target sequence, and record the two matched subsequences as matching anchors; based on the matching anchors, divide the query sequence and the target sequence into intervals respectively, and obtain the alignment interval of the target sequence and the alignment interval of the query sequence; According to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between the alignment intervals, optimize the alignment intervals of the target sequence and the query sequence respectively to obtain the corresponding optimized alignment intervals; Alignment result generation module: Each optimized alignment interval is processed by the double-sequence alignment algorithm to obtain a locally optimal alignment path; according to the positions of the matching anchor points, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

[0013] Advantageous effects: The genomic double-sequence alignment method based on syncmer provided by the present invention first samples based on the syncmer strategy, extracts subsequences, and divides them according to the subsequences, dividing continuous long alignment intervals into multiple short alignment intervals, which significantly improves the computational efficiency while ensuring a relatively high alignment accuracy; and dynamically optimizes and merges the short alignment intervals to effectively balance the interval continuity and alignment efficiency; in order to reduce the influence of interval division on the alignment result, according to the positions of the matching anchor points, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate a more accurate alignment result. Description of the Drawings

[0014] Figure 1 This is a comparison of the running time of a genomic double-sequence alignment method based on syncmer of the present invention with the double-sequence alignment algorithms Ksw2 and Wfa2.

[0015] Figure 2 This is a comparison of the peak memory occupancy of a genomic double-sequence alignment method based on syncmer of the present invention with the double-sequence alignment algorithms Ksw2 and Wfa2.

[0016] Figure 3 This is a comparison of the accuracy of a genomic double-sequence alignment method based on syncmer of the present invention with the double-sequence alignment algorithms Ksw2 and Wfa2. Detailed Embodiments

[0017] The following embodiments are intended to illustrate the present invention rather than further limit the present invention.

[0018] The present invention provides a genomic double-sequence alignment method based on syncmer, including the following steps: S1: Obtain a target sequence and a query sequence, and use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively. If the distance between two adjacent subsequences is greater than the distance threshold, insert at least one subsequence into the interval between the two adjacent subsequences so that the distance between the two adjacent subsequences is not greater than the distance threshold, obtaining the subsequence of the target sequence and the subsequence of the query sequence; After constructing hash indexes based on the subsequence of the target sequence and the subsequence of the query sequence respectively, a target sequence index and a query sequence index are obtained.

[0019] Specifically as follows: First, obtain the sequence information of the target sequence, and extract the sequence name and the base sequence; obtain the sequence information of the query sequence, and extract the sequence name and the base sequence.

[0020] Secondly, regarding syncmer, it is a method used to select conserved k-mers in bioinformatics, and k-mers are selected by examining the positions of substrings of length s (s < k) in the k-mer. Further, use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively, specifically: Select a k-mer within a sliding window containing multiple k-mers (DNA sequences of length k). If the s-mer (subsequence of length s, where s < k) is at the starting position or the end position of the k-mer, then regard the k-mer as a subsequence; continue to slide the current window and extract subsequences until the subsequences of the target sequence and the query sequence are all extracted.

[0021] All in all, extracting the subsequences of the target sequence or the query sequence means extracting a specific (meeting the above conditions) k-mer within each small range, that is, a sliding window.

[0022] Store all the subsequences extracted from the target sequence and their position information in the sequence into the array T[ ]. Store all the subsequences extracted from the query sequence and their position information in the sequence into the array Q[ ].

[0023] Considering that after the initial sampling based on syncmer is completed for the target sequence and the query sequence, there may be a problem that the position distance between adjacent subsequences in the sequence is too large, which may lead to fragmentation of subsequent interval division and thus reduce the calculation efficiency. For this reason, the present invention proposes a dynamic secondary sampling method to optimize the distribution uniformity.

[0024] Specifically, set a spacing threshold E for two adjacent subsequences. If the spacing between two adjacent subsequences is greater than the spacing threshold, then insert at least one subsequence into the interval between the two adjacent subsequences, which can be achieved by relaxing the selection conditions of the subsequences (such as reducing the minimum s-mer strictness or adjusting the hash screening threshold) so that the spacing between the two adjacent subsequences is not greater than the spacing threshold.

[0025] During the secondary sampling process, strictly retain the entire subsequence set of the initial sampling, and only supplement the newly added subsequences, so as to ensure that the subsequence sets of the target sequence and the query sequence are completely compatible after the secondary sampling. This strategy significantly improves the distribution uniformity of subsequences in the sequence by dynamically filling sparse regions, makes the subsequent alignment interval division more coherent, and at the same time avoids redundant calculations caused by sparse distribution, overall improving the algorithm efficiency and result reliability.

[0026] Next, to accelerate the query speed, a hash index was constructed. Specifically, the base sequence of the subsequence of the target sequence was used as the key of the hash table, and the position where the subsequence of the target sequence appears in the target sequence was used as the value of the hash table. The base sequence of the subsequence of the query sequence was used as the key of the hash table, and the position where the subsequence of the query sequence appears in the query sequence was used as the value of the hash table, so as to store the data in the arrays T[] and Q[] into the hash table hash.

[0027] Through the above operations, a genomic index was obtained: the target sequence index and the query sequence index, which include the base sequences of the target sequence and the query sequence, as well as all their subsequence sets T[] and Q[], and the corresponding hash index hash.

[0028] S2: According to the target sequence index, the query sequence index, and each subsequence of the query sequence, match in the target sequence, and mark the two successfully matched subsequences as matching anchors; based on the matching anchors, divide the query sequence and the target sequence into intervals respectively to obtain the alignment interval of the target sequence and the alignment interval of the query sequence; According to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between the alignment intervals, optimize the alignment intervals of the target sequence and the query sequence respectively to obtain the corresponding optimized alignment intervals.

[0029] In this stage, to solve the problem of excessive consumption of time and space resources in the double-sequence alignment process of long reads, the operations are as follows: First, according to the target sequence index, the query sequence index, and each subsequence of the query sequence, match in the target sequence, and mark the two matched subsequences as matching anchors. Specifically: Using the base sequence of the subsequence of the query sequence as the key, search in the hash index of the target sequence. If a matching subsequence is found, mark the two matched subsequences as matching anchors, and record the base sequence of this subsequence, the position of this subsequence in the target sequence, and the position of this subsequence in the query sequence.

[0030] In actual operation, each subsequence is sequentially extracted from the subsequence array Q[] of the query sequence. For each subsequence in the query sequence, search for a match in the hash index of the target sequence: use the base sequence of the subsequence of the query sequence as the key to search in the hash table of the target sequence. If a matching subsequence is found, record the successfully matched subsequence and its position information, including the base sequence of this subsequence, the position (pos_target) of this subsequence in the target sequence, and the position (pos_query) of this subsequence in the query sequence. And store this information into the anchor array C[] as the basis for subsequent interval division.

[0031] Then, based on the matching anchors, the query sequence and the target sequence are respectively divided into intervals to obtain the alignment intervals of the target sequence and the query sequence. Specifically: Sort the matching anchors of the query sequence and the target sequence respectively. Taking two adjacent anchors as boundaries, after dividing the query sequence or the target sequence, use the termination position of the former of the two adjacent anchors as the starting position of the query sequence or the target sequence, and use the starting position of the latter of the two adjacent anchors as the termination position of the query sequence or the target sequence to obtain the alignment intervals of the target sequence and the query sequence.

[0032] However, during the process of interval division, the short interval division effect is often not ideal due to sparse distribution or blurred boundaries of adjacent intervals. For example, it is unable to effectively handle redundant divisions of overlapping intervals and short interval intervals, or wrongly merges some discontinuous regions and other problems. To address this issue, the present invention proposes a method for dynamically optimizing and merging short intervals after subsequence division, which is used to dynamically balance interval continuity and improve alignment efficiency.

[0033] Preferably, according to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between alignment intervals, the alignment intervals of the target sequence and the query sequence are respectively optimized to obtain the corresponding optimized alignment intervals. Specifically: For the target sequence or the query sequence, if the distance between two adjacent alignment intervals is less than the preset threshold D, then the two adjacent alignment intervals are merged to obtain the corresponding optimized alignment interval; wherein, the preset threshold is set according to the sequence characteristics of the alignment interval to be optimized (such as GC offset, base conservation) to avoid loss of alignment information caused by rigid cutting and division; If the length of the alignment interval is less than the preset length (such as less than 100bp), that is, the interval is not significant, then the alignment interval is merged with the adjacent alignment interval to obtain the corresponding optimized alignment interval; If there is an overlap between alignment intervals (such as complete overlap or nesting), then the overlapping alignment intervals are merged and redundant boundaries are removed to obtain the corresponding optimized alignment interval.

[0034] This process significantly improves the continuity and computational efficiency of the alignment by dynamically adjusting the interval boundaries, and avoids performance loss caused by fragmented intervals.

[0035] Furthermore, the matching anchors can be filtered, and low-quality matches are removed according to the matching score (such as position offset or sequence similarity), and the best matching position of each subsequence in the target sequence is retained to ensure the reliability of the alignment interval.

[0036] Furthermore, perform a secondary check on the optimized alignment regions. If the distance between adjacent optimized alignment regions is less than the threshold E, further merge them to obtain the final alignment region list Final_Regions[ ], providing a high-confidence input for subsequent base-level fine-grained alignment.

[0037] S3: Process each optimized alignment region through a pairwise alignment algorithm to obtain a locally optimal alignment path; according to the positions of the matching anchors, splice the locally optimal alignment paths of adjacent optimized alignment regions to generate an alignment result.

[0038] In this stage, to reduce the impact of the S2 partitioning operation on the S3 alignment result, preferably, according to the positions of the matching anchors, splice the locally optimal alignment paths of adjacent optimized alignment regions to generate an alignment result. Specifically: According to the positions of the matching anchors, if the offsets of the locally optimal alignment paths of adjacent optimized alignment regions on the target sequence and the query sequence are the same, merge the locally optimal alignment paths of adjacent optimized alignment regions into a continuous alignment path; If there are overlaps or contradictions in the locally optimal alignment paths of adjacent optimized alignment regions, adopt a dynamic programming backtracking strategy to generate a globally optimal path based on the weights of the locally optimal alignment paths; After splicing the locally optimal alignment paths of adjacent optimized alignment regions corresponding to all matching anchor positions, generate an alignment result.

[0039] In actual operation, use the optimized alignment regions in the S2 alignment region list Final_Regions[ ] as input, and call an efficient pairwise alignment algorithm (such as Ksw2, Wfa2) for parallel processing. Each optimized alignment region independently calculates the locally optimal alignment path and score, and saves the alignment result.

[0040] When splicing the locally optimal alignment paths of adjacent optimized alignment regions according to the position information of the matching anchors in sequence continuity, the scores of each optimized alignment region are weighted and summed according to the coverage length, continuously integrating the locally optimal alignment paths of each anchor and its adjacent optimized alignment regions and merging the optimal alignment scores. Among them, if the offsets of the locally optimal alignment paths of adjacent optimized alignment regions on the target sequence and the query sequence are consistent (pos_target = pos_query), directly merge them; for overlapping or contradictory paths, that is, there are multiple high-score paths in the same region, adopt a dynamic programming backtracking strategy to select the globally optimal path based on the weights of the locally optimal alignment paths.

[0041] In addition, perform base-level fine-grained alignment on sequence fragments that have not been aligned or covered by anchors, and re-integrate the final optimal alignment path and merge the optimal alignment scores to generate a final continuous alignment path.

[0042] Finally, the output alignment results include the global alignment score, the global alignment path, the alignment interval, and the statistical information of the alignment results, such as the number of matching bases, the number of mismatched bases, and the number of inserted or deleted bases.

[0043] The genomic dual-sequence alignment method based on syncmer provided by the present invention first samples based on the syncmer strategy, extracts subsequences, and divides them according to the subsequences, dividing continuous long alignment intervals into multiple short alignment intervals, which significantly improves the computational efficiency while ensuring a high alignment accuracy; and dynamically optimizes and merges the short alignment intervals, effectively balancing the interval continuity and the alignment efficiency; in order to reduce the influence of interval division on the alignment results, according to the positions of the matching anchors, the local optimal alignment paths of adjacent optimized alignment intervals are spliced to generate a more accurate alignment result.

[0044] The present invention also provides a genomic dual-sequence alignment system based on syncmer, including: Index construction module: used to obtain the target sequence and the query sequence, extract the subsequences of the target sequence and the query sequence respectively using the syncmer strategy, and if the distance between two adjacent subsequences is greater than the distance threshold, insert at least one subsequence between the two adjacent subsequence intervals so that the distance between the two adjacent subsequences is not greater than the distance threshold, to obtain the subsequence of the target sequence and the subsequence of the query sequence; After constructing hash indexes respectively based on the subsequence of the target sequence and the subsequence of the query sequence, obtain the target sequence index and the query sequence index; Alignment interval optimization module: according to the target sequence index and the query sequence index, each subsequence of the query sequence is matched in the target sequence, and the two matched subsequences are recorded as matching anchors; based on the matching anchors, the target sequence and the query sequence are respectively divided into intervals to obtain the alignment intervals of the target sequence and the alignment intervals of the query sequence; According to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between the alignment intervals, the alignment intervals of the target sequence and the alignment intervals of the query sequence are respectively optimized to obtain the corresponding optimized alignment intervals; Alignment result generation module: each optimized alignment interval is processed by the dual-sequence alignment algorithm to obtain the local optimal alignment path; according to the positions of the matching anchors, the local optimal alignment paths of adjacent optimized alignment intervals are spliced to generate the alignment result.

[0045] Experiment 1: Comparative experiment on the running speed of the method of the present application and the existing dual-sequence alignment algorithm To verify the improvement in the alignment speed of the method of the present application, the widely used existing pairwise alignment algorithms Ksw2 and Wfa2 were selected for comparison in the experiment. The same simulated dataset was used in the experiment, and an affine gap penalty strategy more representative of real biological sequences was adopted to more accurately simulate the costs of insertions and deletions in biology. In this experiment, the base match score was set to +2, the mismatch score was set to -2, the gap opening score was set to -4, and the gap extension score was set to -2. The experiment was run on a T640 server with an Intel Xeon Gold 6226r processor (main frequency 2.9 GHZ, 16 cores), 384 GB of memory, and an Ubuntu 20.04 operating system. The specific experimental steps are as follows: (1) Use the simulated sequence tool Pbsim3 to simulate sequencing sequences with a depth of 5 times for chromosome 1 of the human reference genome GRch38 with default parameters. By intercepting the sequencing sequences and their reference sequences and then processing them, the pairwise sequence data used in the experiment was generated. To obtain more comprehensive and authoritative comparison results, data with three different lengths (50k, 100k, 200k) and two different error rates (1%, 5%) were generated here.

[0046] (2) Write a program to run the method of the present application (denoted as Synpair), Ksw2, and Wfa2 using the pairwise sequence data generated above as input and output the alignment results, including the alignment score and the alignment path. In addition, the resource utilization information was also recorded simultaneously.

[0047] The alignment running speeds of the three alignment methods for data with different lengths and error rates are shown in Table 1, Figure 1 as follows.

[0048] Table 1 Running speeds of the three alignment methods for data with different lengths and error rates (unit: seconds)

[0049] The experimental results show that the method of the present application (Synpair) has a significantly faster running speed than the existing mainstream algorithms Ksw2 and Wfa2 in scenarios of long read length data with different lengths (50k, 100k, 200k) and different error rates (1%, 5%). Among them, for the data with a length of 50k and an error rate of 1%, the alignment speed of the method of the present application is 7 times faster than that of the existing mainstream algorithm Wfa2; while for the high-complexity data scenario with a length of 200k and an error rate of 5%, the speed advantage of the method of the present application is further expanded to 40 times.

[0050] Particularly crucial is that the performance of the method of this application is minimally affected by sequence length and error rate. For data with different lengths and error rates, the fluctuation range of the running speed of the method of this application is always lower than 15%, reflecting the strong robustness and scalability of this method. This characteristic stems from an efficient computing architecture based on the syncmer strategy, which quickly locates high-confidence alignment anchors by screening specific k-mers (syncmers), reducing the quadratic complexity O(n²) of global dynamic programming to linear complexity O(n). At the same time, it combines multi-threaded parallel processing of independent short-range tasks to avoid the memory and computing bottlenecks of long-read data. In addition, the syncmer strategy naturally filters high-variation or low-quality regions, reducing invalid alignment calculations, enabling it to maintain stable performance in high-error-rate scenarios and providing a breakthrough solution for large-scale analysis of third-generation sequencing data.

[0051] Experiment 2: Comparative experiment on the peak memory occupancy during the runtime of the method of this application and existing double-sequence alignment algorithms To evaluate the peak memory occupancy of the method of this application, the widely used existing double-sequence alignment algorithms Ksw2 and Wfa2 were selected for comparison in the experiment. The experimental configuration, experimental parameter settings, data source are the same as those in Experiment 1.

[0052] The peak memory occupancy of the three alignment methods for data with different lengths and error rates is shown in Table 2, Figure 2 as follows.

[0053] Table 2 Peak memory occupancy (unit: megabytes) of the three alignment methods for data with different lengths and error rates

[0054] The peak memory occupancy of the method of this application (Synpair) for all types of data is significantly lower than that of Ksw2 and Wfa2. Taking the data with a length of 100k and an error rate of 5% as an example, the peak memory of Synpair is only 8 megabytes, while the peak memories of Wfa2 and Ksw2 are as high as 3960 and 66350 megabytes. It can be seen that the memory efficiency of the method of this application has been greatly improved compared with existing double-sequence alignment methods.

[0055] This is because the existing algorithms Ksw2 and Wfa2 need to construct a global dynamic programming matrix for two long sequences (the space complexity reaches the quadratic level); for example, for 100k * 100k data, the number of cells in the dynamic programming matrix reaches 10 10 orders of magnitude. While the method of this application cuts the long sequence into independent short alignment intervals through anchors, making the average length of each alignment interval less than 300, and the number of cells in the dynamic programming matrix that needs to be constructed for each interval is only 10 4This breakthrough improvement in memory efficiency enables the method of this application to process ultra-long read length data in a resource-constrained environment, solving the technical bottleneck of existing algorithms that cannot be expanded due to memory hardware limitations, and providing a solution to larger genome alignment and analysis problems.

[0056] Experiment 3: Comparative experiment on the accuracy of the method of this application and the existing double sequence alignment algorithm In order to evaluate the alignment accuracy of the method of the present application, the experiment selected the widely used dual sequence alignment algorithms Ksw2 and Wfa2 for comparison. The experimental configuration, experimental parameter settings and data sources are the same as those of Experiment 1. In order to verify the alignment accuracy of the three alignment methods, this experiment adopted and implemented the most original and standard dynamic programming algorithm, and used the alignment score results obtained as the alignment data of the benchmark comparison method. Then the ratio of the same sequence to all sequences in the alignment scores of the three alignment methods and the benchmark comparison method is used as the accuracy evaluation index. In other words, the most original and standard method is used, and its results are used as the benchmark. Then the results of the three algorithms (Synpair), Ksw2, and Wfa2 of the present application are compared with the benchmark results to obtain the accuracy index.

[0057] The comparison accuracy of the three comparison methods when running programs with different lengths and error rate data is shown in Table 3. Figure 3 It is worth mentioning that both Ksw2 and Wfa2 methods belong to the existing accurate alignment algorithms, so their accuracy can reach 100%.

[0058] Table 3 Accuracy of three comparison methods under different lengths and error rates

[0059] The experimental results show that the accuracy of the method (Synpair) of this application can achieve ideal results under different data. Even in the highly complex data scenario with a length of 200k and an error rate of 5%, the alignment accuracy of the method of this application can still reach 93%. Since the method of this application adopts a subsequence extraction strategy to reduce the alignment interval, it replaces the traditional base-level alignment, and achieves a significant improvement in time and space efficiency. At the expense of a part of the alignment accuracy, the time and space efficiency can be improved by several to dozens of times, which is completely acceptable in practical applications. Especially when processing large-scale genomic data, its efficiency advantage can significantly reduce the computing resource overhead and provide an efficient solution for long sequence alignment. This method sacrifices limited precision in exchange for nearly a hundred-fold speed increase and exponential optimization of memory efficiency, making it possible for complex analysis of population genomes.

Claims

1. A syncmer-based genome pairwise sequence alignment method, characterized in that: include: S1: Obtain the target sequence and the query sequence, and use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively. If the distance between two adjacent subsequences is greater than the distance threshold, insert at least one subsequence between the two adjacent subsequences so that the distance between the two adjacent subsequences is not greater than the distance threshold, and obtain the subsequence of the target sequence and the subsequence of the query sequence; Based on the subsequence of the target sequence and the subsequence of the query sequence, hash indexes are constructed respectively to obtain the target sequence index and the query sequence index; S2: According to the target sequence index and the query sequence index, each subsequence of the query sequence is matched in the target sequence, and the two matched subsequences are recorded as matching anchor points; based on the matching anchor points, the query sequence and the target sequence are divided into intervals to obtain the comparison interval of the target sequence and the comparison interval of the query sequence; According to the spacing between two adjacent alignment intervals, the length of the alignment intervals, or the overlap between the alignment intervals, the alignment intervals of the target sequence and the alignment intervals of the query sequence are optimized to obtain corresponding optimized alignment intervals; S3: Each optimized alignment interval is processed by a double sequence alignment algorithm to obtain a local optimal alignment path; according to the matching anchor point position, the local optimal alignment paths of adjacent optimized alignment intervals are spliced ​​to generate an alignment result.

2. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: The step S2 optimizes the alignment interval of the target sequence and the alignment interval of the query sequence according to the spacing between two adjacent alignment intervals, the length of the alignment intervals, or the overlap between the alignment intervals, to obtain corresponding optimized alignment intervals, specifically: For the target sequence or query sequence, if the distance between two adjacent comparison intervals is less than a preset threshold, the two adjacent comparison intervals are merged to obtain a corresponding optimized comparison interval; wherein the preset threshold is set according to the sequence characteristics of the comparison interval to be optimized; If the length of the comparison interval is less than the preset length, the comparison interval is merged with the adjacent comparison interval to obtain the corresponding optimized comparison interval; If there is overlap between the comparison intervals, the overlapping comparison intervals are merged to obtain the corresponding optimized comparison intervals.

3. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: S3, according to the matching anchor point position, splices the local optimal comparison paths of adjacent optimized comparison intervals to generate a comparison result, specifically: According to the matching anchor point position, if the local optimal alignment paths of adjacent optimized alignment intervals have the same offset on the target sequence and the query sequence, the local optimal alignment paths of adjacent optimized alignment intervals are merged into one continuous alignment path; If there are overlaps or contradictions in the local optimal comparison paths of adjacent optimization comparison intervals, a dynamic programming backtracking strategy is used to generate a global optimal path based on the weights of the local optimal comparison paths; The local optimal comparison paths of adjacent optimized comparison intervals corresponding to all matching anchor point positions are spliced ​​together to generate a comparison result.

4. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: The step S2 is to divide the query sequence and the target sequence into intervals based on the matching anchor points to obtain the comparison interval of the target sequence and the comparison interval of the query sequence, which is specifically: The matching anchor points of the query sequence and the target sequence are sorted respectively, and the query sequence or the target sequence is divided with two adjacent anchor points as the boundary. The end position of the former of the two adjacent anchor points is used as the start position of the query sequence or the target sequence, and the start position of the latter of the two adjacent anchor points is used as the end position of the query sequence or the target sequence, to obtain the comparison interval of the target sequence and the comparison interval of the query sequence.

5. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: S2, according to the target sequence index and the query sequence index, matches each subsequence of the query sequence in the target sequence, and records the two successfully matched subsequences as matching anchor points, specifically: The base sequence of the subsequence of the query sequence is used as the key to search in the hash index of the target sequence. If a matching subsequence is found, the two successfully matched subsequences are recorded as matching anchor points, and the base sequence of the subsequence, the position of the subsequence in the target sequence, and the position of the subsequence in the query sequence are recorded.

6. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: S1, using syncmer to extract subsequences of the target sequence and the query sequence respectively, is specifically: Select a k-mer in a sliding window containing multiple k-mers. If the s-mer is at the starting position or the end position of the k-mer, the k-mer is regarded as a subsequence. The current window continues to slide to extract subsequences until the subsequences of the target sequence and the query sequence are extracted.

7. The syncmer-based genome pair sequence alignment method according to claim 1, characterized in that: S1, based on the subsequence of the target sequence and the subsequence of the query sequence, respectively constructs a hash index, wherein the base sequence of the subsequence of the target sequence is used as the key of the hash table, and the position where the subsequence of the target sequence appears in the target sequence is used as the value of the hash table; The base sequence of the subsequence of the query sequence is used as the key of the hash table, and the position where the subsequence of the query sequence appears in the query sequence is used as the value of the hash table.

8. The syncmer-based genome pair sequence alignment method according to claim 3, characterized in that: The comparison results include global comparison score, global comparison path, comparison interval, number of matched bases, number of mismatched bases, and number of inserted or deleted bases.

9. A syncmer-based genome binary sequence alignment system, characterized in that: include: Index construction module: used to obtain the target sequence and the query sequence, and use the syncmer strategy to extract subsequences of the target sequence and the query sequence respectively. If the distance between two adjacent subsequences is greater than the distance threshold, at least one subsequence is inserted between the two adjacent subsequences so that the distance between the two adjacent subsequences is not greater than the distance threshold, and the subsequence of the target sequence and the subsequence of the query sequence are obtained; Based on the subsequence of the target sequence and the subsequence of the query sequence, hash indexes are constructed respectively to obtain the target sequence index and the query sequence index; Alignment interval optimization module: According to the target sequence index and the query sequence index, each subsequence of the query sequence is matched in the target sequence, and the two matched subsequences are recorded as matching anchor points; based on the matching anchor points, the query sequence and the target sequence are divided into intervals to obtain the alignment interval of the target sequence and the alignment interval of the query sequence; According to the spacing between two adjacent alignment intervals, the length of the alignment intervals, or the overlap between the alignment intervals, the alignment intervals of the target sequence and the alignment intervals of the query sequence are optimized to obtain corresponding optimized alignment intervals; Alignment result generation module: Each optimized alignment interval is processed by a double sequence alignment algorithm to obtain a local optimal alignment path; according to the matching anchor point position, the local optimal alignment paths of adjacent optimized alignment intervals are spliced ​​to generate an alignment result.

Images