Syncmer-based genomic dual-sequence alignment method and system

Abstract

The present invention belongs to the technical field of gene sequence alignment, and particularly relates to a genomic dual-sequence alignment method and system based on syncmer. First, sampling is performed based on the syncmer strategy, subsequences are extracted, and division is carried out according to the subsequences, dividing continuous long alignment intervals into multiple short alignment intervals, which significantly improves the calculation efficiency while ensuring relatively high alignment accuracy; moreover, dynamic optimization and merging are performed on the short alignment intervals, effectively balancing interval continuity and alignment efficiency; in order to reduce the influence of interval division on the alignment result, according to the positions of matching anchor points, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate a more accurate alignment result.

Description

Technical Field

[0001] The present invention belongs to the technical field of gene sequence alignment, and particularly relates to a genomic dual-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, dual-sequence alignment, as the most basic 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 innovation in efficient algorithms. Summary of the Invention

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

[0004] The technical solution of the present invention is as follows:

[0005] The present invention provides a genomic dual-sequence alignment method based on syncmer, including:

[0006] S1: Obtain a target sequence and a query sequence, and respectively extract subsequences of the target sequence and the query sequence by using the syncmer strategy. 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;

[0007] After respectively constructing hash indexes 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;

[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 matched subsequences as matching anchors; respectively divide the intervals of the query sequence and the target sequence based on the matching anchors to obtain the alignment interval of the target sequence and the alignment interval of the query sequence;

[0009] According to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between the alignment intervals, respectively optimize the alignment intervals of the target sequence and the alignment intervals of the query sequence to obtain the corresponding optimized alignment intervals;

[0010] S3: Each optimized alignment interval is processed by a pairwise alignment algorithm to obtain a locally optimal alignment path; according to the positions of the matching anchors, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

[0011] In step S2, 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 query sequence are optimized respectively to obtain corresponding optimized alignment intervals. Specifically:

[0012] For the target sequence or the query sequence, if the distance between two adjacent alignment intervals is less than a preset threshold, the two adjacent alignment intervals are merged to obtain a corresponding optimized alignment interval; where the preset threshold is set according to the sequence characteristics of the alignment interval to be optimized.

[0013] If the length of the alignment interval is less than a preset length, the alignment interval is merged with an adjacent alignment interval to obtain a corresponding optimized alignment interval.

[0014] If there is an overlap between the alignment intervals, the overlapping alignment intervals are merged to obtain a corresponding optimized alignment interval.

[0015] In step S3, according to the positions of the matching anchors, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result. Specifically:

[0016] According to the positions of the matching anchors, 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, the locally optimal alignment paths of the adjacent optimized alignment intervals are merged into a continuous alignment path.

[0017] If there is an overlap or contradiction in the locally optimal alignment paths of adjacent optimized alignment intervals, a dynamic programming backtracking strategy is adopted to generate a globally optimal path based on the weights of the locally optimal alignment paths.

[0018] After splicing the locally optimal alignment paths of adjacent optimized alignment intervals corresponding to all the positions of the matching anchors, an alignment result is generated.

[0019] In step S2, the query sequence and the target sequence are respectively partitioned into intervals based on the matching anchors to obtain the alignment intervals of the target sequence and the query sequence. Specifically:

[0020] The matching anchors of the query sequence and the target sequence are sorted respectively. Taking two adjacent anchors as boundaries, after partitioning the query sequence or the target sequence, the start position of the query sequence or the target sequence is the termination position of the former of the two adjacent anchors, and the termination position of the query sequence or the target sequence is the start position of the latter of the two adjacent anchors, so as to obtain the alignment intervals of the target sequence and the query sequence.

[0021] S2 queries each subsequence of the query sequence in the target sequence according to the target sequence index and the query sequence index, and records the two successfully matched subsequences as matching anchors. Specifically:

[0022] 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 this subsequence, the position of this subsequence in the target sequence, and the position of this subsequence in the query sequence.

[0023] S1 uses syncmers to extract subsequences of the target sequence and the query sequence respectively. Specifically:

[0024] 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, regard 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 all extracted.

[0025] S1 constructs hash indexes respectively based on the subsequences of the target sequence and the subsequences of the query sequence. Use the base sequence of the subsequence of the target sequence as the key of the hash table, and the position where the subsequence of the target sequence appears in the target sequence as the value of the hash table;

[0026] Use the base sequence of the subsequence of the query sequence as the key of the hash table, and the position where the subsequence of the query sequence appears in the query sequence as the value of the hash table.

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

[0028] The present invention also provides a genomic double-sequence alignment system based on syncmers, including:

[0029] 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, to obtain the subsequences of the target sequence and the subsequences of the query sequence;

[0030] After constructing hash indexes respectively based on the subsequences of the target sequence and the subsequences of the query sequence, obtain the target sequence index and the query sequence index;

[0031] Alignment interval optimization module: According to the target sequence index and the query sequence index, each subsequence of the query sequence is searched and matched in the target sequence, and the two matched subsequences are recorded as matching anchors; based on the matching anchors, the query sequence and the target sequence are respectively divided into intervals to obtain the alignment interval of the target sequence and the alignment interval of the query sequence;

[0032] 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;

[0033] 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 position of the matching anchor, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

[0034] Beneficial 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 calculation efficiency while ensuring a high alignment accuracy; and dynamically optimizes and merges the short alignment intervals to effectively balance the interval continuity and the alignment efficiency; in order to reduce the impact of interval division on the alignment result, according to the position of the matching anchor, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate a more accurate alignment result. Description of the Drawings

[0035] 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.

[0036] 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.

[0037] 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

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

[0039] The present invention provides a genomic double-sequence alignment method based on syncmer, including the following steps:

[0040] 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 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 subsequences of the target sequence and the subsequences of the query sequence.

[0041] After constructing hash indexes based on the subsequences of the target sequence and the subsequences of the query sequence respectively, obtain the target sequence index and the query sequence index.

[0042] Specifically as follows:

[0043] 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.

[0044] Secondly, regarding syncmer, it is a method used to select conserved k-mers in bioinformatics, and k-mers are selected by checking 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:

[0045] 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 extracted completely.

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

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

[0048] 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 positions of adjacent subsequences in the sequence are too far apart, 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.

[0049] Specifically, a spacing threshold E for two adjacent subsequences is set. If the spacing between two adjacent subsequences is greater than the spacing threshold, at least one subsequence is inserted 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.

[0050] During the secondary sampling process, all subsequence sets of the primary sampling are strictly retained, and only newly added subsequences are supplemented, thus ensuring 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, making the subsequent alignment interval division more coherent, and at the same time avoiding redundant calculations caused by sparse distribution, thereby overall improving the algorithm efficiency and result reliability.

[0051] Next, to accelerate the query speed, a hash index is constructed. Specifically, 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, so as to store the data in the arrays T[] and Q[] into the hash table hash.

[0052] Through the above operations, a genomic index is 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.

[0053] 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 successfully matched subsequences are recorded as matching anchors; based on the matching anchors, the query sequence and the target sequence are respectively divided into intervals to obtain the alignment interval of the target sequence and the alignment interval of the query sequence;

[0054] According to the spacing 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.

[0055] 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:

[0056] First, 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. Specifically:

[0057] Using the base sequence of a subsequence of the query sequence as a key, search in the hash index of the target sequence. If a matching subsequence is found, record the two matching 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.

[0058] In actual operation, extract each subsequence from the subsequence array Q[] of the query sequence in turn. 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 a key and 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 in the anchor array C[] as the basis for subsequent interval division.

[0059] Then, based on the matching anchors, divide the query sequence and the target sequence into intervals respectively to obtain the alignment intervals of the target sequence and the alignment intervals of the query sequence. Specifically:

[0060] 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 alignment intervals of the query sequence.

[0061] 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 solve this problem, 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.

[0062] Preferably, optimize the alignment intervals of the target sequence and the alignment intervals 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. Specifically:

[0063] For the target sequence or the query sequence, if the distance between two adjacent alignment intervals is less than the preset threshold D, 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 (such as GC offset, base conservation) of the alignment interval to be optimized to avoid loss of alignment information caused by rigid cutting and division;

[0064] If the length of the alignment interval is less than the preset length (e.g., 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;

[0065] If there is an overlap between the 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.

[0066] 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.

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

[0068] Furthermore, a secondary check is performed on the optimized alignment interval. If the distance between adjacent optimized alignment intervals is less than the threshold E, then further merging is performed to obtain the final alignment interval list Final_Regions[ ], providing a high-confidence input for the subsequent base-level fine alignment.

[0069] S3: Each optimized alignment interval is processed by a double-sequence alignment algorithm to obtain a locally optimal alignment path; according to the matching anchor point positions, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

[0070] In this stage, in order to reduce the impact of the S2 partitioning operation on the S3 alignment result, preferably, according to the matching anchor point positions, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result. Specifically:

[0071] According to the matching anchor point positions, 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 the locally optimal alignment paths of adjacent optimized alignment intervals are merged into a continuous alignment path;

[0072] If there is an overlap or contradiction in the locally optimal alignment paths of adjacent optimized alignment intervals, then a dynamic programming backtracking strategy is adopted to generate a globally optimal path based on the weights of the locally optimal alignment paths;

[0073] After splicing the locally optimal alignment paths of adjacent optimized alignment intervals corresponding to all matching anchor point positions, an alignment result is generated.

[0074] In actual operation, the optimized alignment intervals in the S2 comparison interval list Final_Regions[] are used as input, and efficient pairwise alignment algorithms (such as Ksw2, Wfa2) are called for parallel processing. Each optimized alignment interval independently calculates the locally optimal alignment path and score, and saves the alignment result.

[0075] Using the position information of the matching anchors, when splicing the locally optimal alignment paths of adjacent optimized alignment intervals according to sequence continuity, the scores of each optimized alignment interval are weighted and summed according to the coverage length, continuously integrating the locally optimal alignment paths of each anchor and its adjacent optimized alignment intervals and combining the optimal alignment scores. Among them, 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 (pos_target = pos_query), they are directly merged; for overlapping or conflicting paths, that is, there are multiple high-score paths in the same region, a dynamic programming backtracking strategy is adopted to select the globally optimal path based on the weights of the locally optimal alignment paths.

[0076] In addition, for sequence fragments that have not been aligned or are not covered by anchors, a base-level fine alignment is performed, and the final optimal alignment path and the combined optimal alignment score are integrated again to generate the final continuous alignment path.

[0077] 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.

[0078] The syncmer-based pairwise genome alignment method 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 impact of interval division on the alignment result, according to the position of the matching anchors, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate a more accurate alignment result.

[0079] The present invention also provides a syncmer-based pairwise genome alignment system, including:

[0080] 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, if the distance between two adjacent subsequences is greater than the distance threshold, at least one subsequence is inserted 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;

[0081] After constructing hash indexes based on subsequences of the target sequence and subsequences of the query sequence respectively, a target sequence index and a query sequence index are obtained.

[0082] 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 query sequence and the target sequence are respectively divided into intervals to obtain the alignment interval of the target sequence and the alignment interval of the query sequence.

[0083] According to the spacing 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 alignment intervals of the query sequence are respectively optimized to obtain corresponding optimized alignment intervals.

[0084] Alignment result generation module: 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, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

[0085] Experiment 1: A comparison experiment on the running speed of the method of this application and existing double-sequence alignment algorithms

[0086] To verify the improvement in the alignment speed of the method of this application, the existing widely used double-sequence alignment algorithms Ksw2 and Wfa2 were selected for comparison in the experiment. The experiment used the same simulated data set and adopted an affine gap penalty strategy that is more representative of real biological sequences 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.9GHZ, 16 cores), 384GB of memory, and an Ubuntu 20.04 operating system. The specific experimental steps are as follows:

[0087] (1) Use the simulation 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 double-sequence data used in the experiment was generated. To obtain more comprehensive and authoritative comparison results, data of three different lengths (50k, 100k, 200k) and two different error rates (1%, 5%) were generated here.

[0088] (2)Write a program to make the method of this application (denoted as Synpair) and Ksw2, Wfa2, run the generated double-sequence data as input and output the alignment results, including alignment scores and alignment paths. In addition, resource utilization information is also recorded simultaneously.

[0089] The alignment running speeds of the three alignment methods under different lengths and error rate data are shown in Table 1, Figure 1 as follows.

[0090] Table 1 Running speeds of the three alignment methods under different lengths and error rate data (unit: seconds)

[0091]

[0092] The experimental results show that the method of this application (Synpair) has significantly better running speeds than the existing mainstream algorithms Ksw2 and Wfa2 in scenarios of long-read 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 this 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 this application is further expanded to 40 times.

[0093] 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 the 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 the linear level 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.

[0094] Experiment 2: Comparative experiment on the peak memory occupancy during the running of the method of this application and existing double-sequence alignment algorithms

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

[0096] The peak memory usage of the three alignment methods under different lengths and error rate data is shown in Table 2, Figure 2 as shown below.

[0097] Table 2 Peak memory usage (unit: megabytes) of the three alignment methods under different lengths and error rate data

[0098]

[0099] The peak memory usage of the method (Synpair) of this application is significantly lower than that of Ksw2 and Wfa2 under various data. 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 the existing pairwise alignment methods.

[0100] 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 data of 100k * 100k, the number of cells in the dynamic programming matrix reaches 10 10 orders of magnitude. However, the method of this application cuts the long sequences into independent short alignment intervals through anchor points, so that the average length of each alignment interval is less than 300, and the number of cells in the dynamic programming matrix to be constructed for each interval is only 10 4 orders of magnitude. This breakthrough improvement in memory efficiency enables the method of this application to process ultra-long read data in resource-constrained environments, solves the technical bottleneck that the existing algorithms cannot be extended due to memory hardware limitations, and provides a solution to larger genome alignment analysis problems.

[0101] Experiment 3: Accuracy comparison experiment between the method of this application and existing pairwise alignment algorithms

[0102] To evaluate the alignment accuracy of the method of this application, the existing widely used pairwise alignment algorithms Ksw2 and Wfa2 were selected for comparison in the experiment. The experimental configuration, experimental parameter settings, data sources are the same as those in Experiment 1. To verify the alignment accuracy of the three alignment methods, the most primitive and standard dynamic programming algorithm was adopted and implemented in this experiment, and the alignment score results obtained by it were used as the alignment data of the benchmark comparison method. Then, the ratio of the same sequences to all sequences in the alignment scores of the three alignment methods and the benchmark comparison method was used as the accuracy evaluation index. In other words, the most primitive and standard method was used, and its result was used as the benchmark. Then, the results of the three algorithms of the method of this application (Synpair), Ksw2, and Wfa2 were compared with the benchmark result to obtain the accuracy index.

[0103] 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%.

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

[0105]

[0106] 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 genomic dual-sequence alignment method, characterized in that, 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 in 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, 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 anchors; based on the matching anchors, the query sequence and the target sequence are respectively divided into intervals, obtaining 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 overlap 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; 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, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced to generate an alignment result.

2. The syncmer-based genomic dual sequence alignment method according to claim 1, wherein In the said S2, according to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlap 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, specifically: For the target sequence or the query sequence, if the distance between two adjacent alignment intervals is less than the preset threshold, 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; If the length of the alignment interval is less than the preset length, the alignment interval is merged with the adjacent alignment interval to obtain the corresponding optimized alignment interval; If there is an overlap between the alignment intervals, the overlapping alignment intervals are merged to obtain the corresponding optimized alignment interval.

3. The syncmer-based genomic double sequence alignment method according to claim 1, wherein In the said S3, according to the positions of the matching anchors, the locally optimal alignment paths of adjacent optimized alignment intervals are spliced 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 intervals on the target sequence and the query sequence are the same, the locally optimal alignment paths of the adjacent optimized alignment intervals are merged into a continuous alignment path; If there is an overlap or contradiction in the locally optimal alignment paths of adjacent optimized alignment intervals, a dynamic programming backtracking strategy is adopted 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 the positions of the matching anchors, an alignment result is generated.

4. The syncmer-based genomic dual sequence alignment method according to claim 1, wherein In the said S2, based on the matching anchors, the query sequence and the target sequence are respectively divided into intervals to obtain the alignment interval of the target sequence and the alignment interval of the query sequence, specifically: Sort the matching anchors of the query sequence and the target sequence respectively. After dividing the query sequence or the target sequence with two adjacent anchors as the boundaries, 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 interval of the target sequence and the alignment interval of the query sequence.

5. The syncmer-based genomic double sequence alignment method according to claim 1, wherein In step 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 successfully matched subsequences are recorded as matching anchors. Specifically: Use the base sequence of the subsequence of the query sequence as the key to 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 this subsequence, the position of this subsequence in the target sequence, and the position of this subsequence in the query sequence.

6. The syncmer-based genomic dual sequence alignment method according to claim 1, wherein In step S1, use syncmer 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, regard 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 extracted completely.

7. The syncmer-based genomic dual sequence alignment method according to claim 1, characterized in that In step S1, based on the subsequences of the target sequence and the subsequences of the query sequence, construct hash indexes respectively. Use the base sequence of the subsequence of the target sequence as the key of the hash table, and use the position where the subsequence of the target sequence appears in the target sequence as the value of the hash table; Use the base sequence of the subsequence of the query sequence as the key of the hash table, and use the position where the subsequence of the query sequence appears in the query sequence as the value of the hash table.

8. The syncmer-based genomic dual sequence alignment method according to claim 3, wherein 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.

9. A syncmer-based genomic double sequence alignment system, characterized in that, It includes: 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 to obtain the subsequences of the target sequence and the subsequences of the query sequence; After constructing hash indexes based on the subsequences of the target sequence and the subsequences of 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, 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, 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; Optimize the alignment intervals of the target sequence and the query sequence respectively according to the distance between two adjacent alignment intervals, the length of the alignment interval, or the overlapping situation between the alignment intervals 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.

Images