Multi-sequence parallel comparison method based on maximum precision matching

Abstract

The invention discloses a multi-sequence parallel comparison method based on maximum precision matching, and relates to the technical field of crossing of bioinformatics and computational biology. The method comprises the following steps: splicing character strings of a plurality of sequences to obtain a long sequence, and determining a suffix array and a longest common prefix array of the long sequence by adopting a mode of combining a bucket dividing strategy and induced sequencing; according to the suffix array, searching a maximum accurate matching fragment from the longest common prefix array, and screening an ordered and non-overlapped maximum accurate matching fragment set through a dynamic programming technology; expanding each maximum accurate matching fragment in the maximum accurate matching fragment set to all sequences to form a chain; segmenting the long sequence through a chain to obtain a plurality of sub-sequence sets; and calling third-party multi-sequence comparison software to compare the plurality of sub-sequence sets, and splicing the comparison result of each sub-sequence set to obtain a multi-sequence comparison result. According to the method, the time complexity of ultra-long or large-scale sequence data comparison is reduced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bioinformatics and computational biology, and particularly relates to a multiple sequence parallel alignment method based on maximum exact match. BACKGROUND

[0002] In the evolution process of organisms, some important genes are retained, such as some genes for controlling structure and function, and some genes are mutated. Therefore, if two groups of gene sequences are highly similar, whether the two groups of sequences have a common ancestor can be analyzed; in addition, if a similar sequence to an unknown sequence is found in a known sequence set, the structure or function of the unknown sequence can be predicted. Therefore, the similarity size analysis of multiple sequences is the first step of many biological analyses, and more abundant biological information can be obtained according to the differences and connections between gene sequences, thereby providing an important basis for understanding the evolutionary relationship of organisms, the structure and function of biological macromolecules.

[0003] Multiple sequence alignment (MSA) as a method that can directly find the similarity size of multiple biological sequences has always been one of the research hotspots in the field of bioinformatics. Multiple sequence alignment aligns three or more sequences, and aligns the common parts between the sequences. These aligned parts may be the core gene parts that play a key role in the evolution process.

[0004] However, the traditional multiple sequence alignment algorithm usually has a high computational complexity, and it is often difficult to obtain efficient alignment results when dealing with ultra-long or large-scale sequence data. SUMMARY

[0005] Therefore, it is necessary to provide a multiple sequence parallel alignment method based on maximum exact match, which reduces the time complexity of ultra-long or large-scale sequence data alignment and improves the alignment efficiency.

[0006] The present application adopts the following technical solutions: The present application provides a multiple sequence parallel alignment method based on maximum exact match, comprising: The strings of multiple sequences are spliced to obtain long sequences, and a bucketing strategy and an induced sorting combination are used to determine the suffix array and the longest common prefix array of the long sequences; According to the suffix array, maximum exact match segments are found from the longest common prefix array, and a dynamic programming technique is used to screen an ordered and non-overlapping maximum exact match segment set; Each maximum exact match segment in the maximum exact match segment set is extended to all sequences to form a chain; The long sequence is cut into multiple subsequence sets by chain; The third-party multi-sequence alignment software is called to align the multiple subsequence sets, and the alignment results of each subsequence set are spliced to obtain a multi-sequence alignment result.

[0007] Optionally, there is a separator “$” between different sequences in the long sequence, and there is an end symbol “ ” at the end of the long sequence. In the dictionary order, it is defined as .

[0008] Optionally, a combination of the bucketing strategy and the induced sorting is used to determine the suffix array and the longest common prefix array of the long sequence, including: According to the frequency of characters in the long sequence, the initial buckets are divided; The suffix sorting of special characters is recursively processed, and the suffix array, the longest common prefix array and the document array are improved through two rounds of induced sorting.

[0009] Optionally, the suffix array includes the suffix start position; according to the suffix array, the maximum exact match segment is found from the longest common prefix array, including: The longest common prefix array is traversed to determine the longest common prefix interval; An initialized position list is obtained, and each longest common prefix interval is traversed to add the suffix start position of each longest common prefix interval to the position list; Through left extension, it is checked whether the characters after the suffix start position minus the offset are consistent, and if so, the suffix start position of the longest common prefix interval is extended to the left by offset until it cannot be extended; According to the extended start position and the length of the common prefix, the maximum exact match segment is determined.

[0010] Optionally, the method further includes: If the length of the maximum exact match segment is greater than the preset sequence length, the maximum exact match segment is removed.

[0011] Optionally, the ordered and non-overlapping maximum exact match segment set is screened through dynamic programming technology, including: For any maximum exact match segment, the average position of each substring in the maximum exact match segment is calculated; If there are repeated substrings in the same maximum exact match segment, the substring farthest from the average position is deleted; All maximum exact match segments are arranged in ascending order according to the average position; The conflict problem of all maximum exact match segments is optimized through dynamic programming strategy to obtain an ordered and non-overlapping maximum exact match segment set.

[0012] Optionally, the dynamic programming strategy comprises a fast mode and an accurate mode; the fast mode comprises taking the maximum accurate matching area as a target function, and screening a maximum accurate matching set maximizing the total area by one global dynamic programming; the accurate mode comprises independently running dynamic programming on each maximum accurate matching, taking the maximum fragment length as a target function, and eliminating overlaps in sequence order to retain the maximum accurate matching satisfying the minimum number of covered sequences. M

[0013] Optionally, each maximum accurate matching fragment in the maximum accurate matching fragment set is extended to all sequences to form a chain, comprising: Each maximum accurate matching fragment in the maximum accurate matching fragment set is extended to all sequences to form a chain by an SSW algorithm.

[0014] Optionally, third-party multi-sequence alignment software is called to align the multiple sub-sequence sets, comprising: The G-INS-i algorithm of MAFFT and the fast mode of HAlign3 are used to align the multiple sub-sequence sets.

[0015] Optionally, the method further comprises: A Profile alignment algorithm based on fast Fourier transform and K-Band constraint is used to align the maximum accurate matching fragment that cannot be extended; The alignment results of each sub-sequence set and the alignment results of the maximum accurate matching fragment that cannot be extended are spliced to obtain a multi-sequence alignment result.

[0016] The application provides a multi-sequence parallel alignment device based on maximum accurate matching, comprising: A construction module is configured to splice character strings of multiple sequences to obtain a long sequence, and determine a suffix array and a longest common prefix array of the long sequence by combining a bucketing strategy and an induced sorting method; A screening module is configured to find maximum accurate matching fragments from the longest common prefix array according to the suffix array, and screen an ordered and non-overlapping maximum accurate matching fragment set by a dynamic programming technique; An extension module is configured to extend each maximum accurate matching fragment in the maximum accurate matching fragment set to all sequences to form a chain; An alignment module is configured to split the long sequence by the chain to obtain multiple sub-sequence sets, call third-party multi-sequence alignment software to align the multiple sub-sequence sets, and splice the alignment results of each sub-sequence set to obtain a multi-sequence alignment result.

[0017] ​The application provides a computer readable storage medium, which stores a computer program, and the computer program realizes the maximum exact match based multi-sequence parallel alignment method when executed by a processor.

[0018] The application provides a computer device, which comprises a memory, a processor and a computer program stored in the memory and executable on the processor, and the processor realizes the maximum exact match based multi-sequence parallel alignment method when executing the program.

[0019] The application adopts the above at least one technical scheme to achieve the following beneficial effects: The application combines the bucketing strategy with the induced sorting to construct the suffix array and the longest common prefix array of the long sequence. The bucketing strategy can divide the super-long sequence into multiple local processing sections (buckets) according to the character frequency, thereby reducing the calculation and memory overhead of one-time sorting. The induced sorting effectively guides the relative order calculation of the suffixes on this basis, avoids repeated traversal of the complete suffixes, and thereby realizes the linear time complexity O( N ) of the suffix array construction. N After the construction is completed, the application extracts the maximum exact match fragments, combines the dynamic programming strategy to screen the ordered and non-overlapping maximum exact match fragment set, thereby compresses the original full sequence alignment task into the local alignment of the core region, greatly reduces the data processing scale and the alignment calculation amount. In addition, the chain constructed based on the maximum exact match fragments divides the long sequence into multiple sub-sequence sets, and processes the sub-sequence sets through the parallel multi-sequence alignment tool, effectively converts the original serially executed global alignment task into a sub-task that can be processed concurrently, and significantly improves the overall alignment efficiency. Therefore, when facing super-long or large-scale sequence data, the application significantly reduces the time complexity of the alignment process, and improves the overall computing performance and resource utilization. BRIEF DESCRIPTION OF DRAWINGS

[0020] The accompanying drawings, which are included to provide a further understanding of the application, constitute a part of the application and serve to explain the illustrative embodiments of the application and their descriptions, and do not constitute an improper limitation on the application. In the drawings:

[0021] Figure 1 A maximum exact match based multi-sequence parallel alignment method flowchart provided by the application; Figure 2 An extension process diagram of a SWW algorithm provided by the application; Figure 3 A process diagram of K-band in dynamic programming provided by the application; Figure 4Another maximum exact match based multi-sequence parallel alignment method flow chart provided by the application; Figure 5 The computer device for implementing the maximum exact match based multi-sequence parallel alignment method provided by the application. DETAILED DESCRIPTION

[0022] To make the objectives, technical solutions, and advantages of the application clearer, the technical solutions of the application will be described below in conjunction with specific embodiments of the application and corresponding drawings. Obviously, the described embodiments are only some of the embodiments of the application, rather than all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative work fall within the protection scope of the application.

[0023] At present, the existing vertical division method relies on a global chain division strategy, which is prone to failure when processing low-similarity datasets, and fails to fully utilize the advantages of modern parallel computing. In addition, the current maximum exact match (MEMs) screening method lacks an optimization strategy for dynamic programming, resulting in redundant matching and conflicts during the alignment process, further affecting the accuracy and precision of the alignment.

[0024] To overcome these problems, the application provides a maximum exact match based multi-sequence parallel alignment method, which optimizes MEMs screening, local alignment expansion, and parallel computing technology through dynamic programming, not only effectively reducing the occurrence of redundant matching, but also greatly improving the alignment efficiency, solving the alignment problems of large-scale, low-similarity, and ultra-long sequences. This method has significant advantages when processing complex genomic data, can provide more accurate and efficient alignment results, and promotes the progress of genomics research and related fields.

[0025] The technical solutions provided by the embodiments of the application will be described in detail below with reference to the drawings.

[0026] Figure 1 The flow chart of the maximum exact match based multi-sequence parallel alignment method provided by the application specifically includes the following steps: S101, concatenate the strings of multiple sequences to obtain a long sequence, and determine the suffix array and the longest common prefix array of the long sequence by combining the bucketing strategy and the induced sorting.

[0027] There is a separator "$" between different sequences in the long sequence, and there is an end symbol "$ " at the end of the long sequence. Other characters in the long sequence are defined as in lexicographic order.

[0028] Take DNA sequence string as an example, given strings , the characters in the string belong to the character set , the length of the string is , wherein represents the size of the character set, and for DNA sequences, the character set is , and the size of the character set is . These characters are spliced into a string , wherein $ represents the separator between different sequences, is the end identifier of the spliced sequence, and for other characters is defined in lexicographic order as . The total length of the spliced long sequence is .

[0029] On the basis of this splicing, the invention adopts the gsaca-k algorithm (based on the improved induced sorting framework) to construct the suffix array (SA), the longest common prefix (LCP) array and the document array (DA). The suffix array SA is defined as the result of strictly sorting all suffixes in lexicographic order, wherein SA[ i ] represents the starting position of the i th smallest suffix in the string, and needs to satisfy the strict prefix sorting relationship for all suffixes. For example, if SA[ i ]=5 and SA[ j ]=3, it indicates that the suffix starting at position 5 is strictly greater than the suffix starting at position 3 in lexicographic order. The LCP array (Longest Common Prefix array) is generated by comparing adjacent suffixes, that is, LCP[ i ] is the maximum value of the common prefix length of the two suffixes corresponding to SA[ i -1] and SA[i]. The document array DA further records the original sequence belonging of each suffix, that is, the suffix at SA[ i ] position is derived from the k th sequence in the original sequence set (such as a segment of the original sequence), which can be used for maximum exact matching (MEM) search. The core algorithm of constructing the data structure, gsaca-k, optimizes the calculation efficiency through the bucketing strategy and the induced sorting technique, and controls the time complexity to , and the space complexity is compressed to . The specific implementation includes three stages: first, the initial buckets are divided according to the character frequency, then the suffix sorting of special characters (such as separators and end symbols) is recursively processed, and finally the SA, LCP array and DA are gradually improved through two rounds of induced sorting.

[0030] Therefore, in one embodiment, a combination of bucketing strategy and induced sorting is used to determine the suffix array and the longest common prefix array of a long sequence, including: dividing the long sequence into initial buckets based on the character frequency; recursively processing the suffix sorting of special characters; and improving the suffix array, the longest common prefix array, and the document array through two rounds of induced sorting.

[0031] Furthermore, the generated suffix structure-related arrays (SA, LCP arrays, and DA) are stored as binary files, and block compression is used to reduce memory usage, which facilitates subsequent fast reading and traversal.

[0032] S102: Based on the suffix array, find the maximum exact matching segment from the longest common prefix array, and use dynamic programming to filter the ordered and non-overlapping set of the maximum exact matching segments.

[0033] In one embodiment, the suffix array includes the suffix start position; based on the suffix array, finding the maximum exact matching segment from the longest common prefix array includes: traversing the longest common prefix array to determine the longest common prefix interval; obtaining an initialized position list, traversing each longest common prefix interval, and adding the suffix start position of each longest common prefix interval to the position list; checking whether the characters after subtracting the offset from all suffix start positions are consistent by left expansion; if consistent, expanding the suffix start position of the longest common prefix interval to the left by offset until it cannot be expanded further; determining the maximum exact matching segment based on the expanded start position and the length of the common prefix.

[0034] Specifically, find the LCP interval: After obtaining the LCP array, the LCP intervals can be obtained by traversing the LCP array. The suffix array is sorted lexicographically. Adjacent suffixes in the sorted array share a common prefix, and multiple intervals with a common prefix exist. Each LCP value represents the length of the longest common prefix between the current suffix and the previous suffix. Therefore, the intervals with a common prefix mentioned above can be called LCP intervals (LCP-interval), abbreviated as LCP. range or This interval meets the following four requirements:

[0035] (1) (2) (3) (4) in, The longest common prefix (LCP) value for the interval is given. For the high side, It is the lower side.

[0036] Therefore, based on the above four requirements, it can be concluded that for all indexes... The suffix must have a length of The same prefix. Given an LCP array, this invention can solve for all minimum length thresholds. l and minimum number of spanning sequences c The LCP intervals are determined by the LCP array. All LCP intervals that meet the conditions can be obtained from the LCP array, and these LCP intervals are the basis for solving MEMS problems.

[0037] MEMs must satisfy the condition that they are not expandable to the left or right. The common prefix found through the LCP interval only ensures that it is not expandable to the right, but it cannot guarantee whether there are identical strings to its left. Therefore, the LCP interval needs to be expanded to the left to obtain an exact matching fragment that conforms to the MEMs definition. Given the LCP interval, the suffix array SA, the document array DA, and the concatenated string... This process is implemented through an algorithm: First, initialize the position list, traverse each LCP interval, add the corresponding suffix start position to the position, then check whether the characters after subtracting the offset from all positions are consistent through left expansion. If they are consistent, increase the offset until it can no longer be expanded. Finally, add the start position and length of the expanded substring to the MEMs set.

[0038] In one embodiment, if the length of the maximum exact match segment exceeds a preset sequence length, that maximum exact match segment is discarded. It is worth noting that the MEMs obtained above are not exhaustive, as some LCP intervals may contain even longer LCP intervals. These longer intervals are discarded to reduce computational load while still meeting the sequence partitioning requirements. Through this optimization, the computational complexity of the MEMs is significantly reduced while ensuring the sequence partitioning objective is achieved.

[0039] In one embodiment, the method of filtering an ordered and non-overlapping set of maximum exact matches using dynamic programming includes: for any maximum exact match, calculating the average position of each substring in the maximum exact match; if a duplicate substring appears in the same maximum exact match, deleting the substring that deviates the furthest from the average position; arranging all maximum exact match segments in ascending order of average position; and optimizing the conflict problem of all maximum exact match segments using a dynamic programming strategy to obtain an ordered and non-overlapping set of maximum exact matches.

[0040] Specifically, after generating the maximum exact matching segment, a two-step filtering process is required to meet the requirement of ordered and conflict-free alignment. The first step is based on sequence coverage deduplication: for each MEM's internal substrings, the average position of each substring is calculated. If a duplicate substring appears in the same sequence, the substring furthest from the average position is deleted. Then, all MEMs are sorted in ascending order of their average positions to ensure global order. The second step introduces a dynamic programming strategy to optimize the conflict problem, which is divided into two modes: fast mode and exact mode.

[0041] The fast mode involves using the maximum exact matching area as the objective function and employing a global dynamic programming approach to select the conflict-free maximum exact matching set that maximizes the total area. Specifically, in the fast mode, the area of ​​the MEM (i.e., "length × number of covering sequences") is used as the objective function, and a global dynamic programming approach is performed to select the conflict-free MEM set with the largest total area. Let there be a total... k Given a number of MEM segments, already sorted by their starting positions, let dp[ i ] indicates the first i The optimal total area ending with MEM is derived from the following recurrence relation:

[0042] dp[ i ] = max(dp[ j ] + len[ i ] × cov[ i ]); in, j < i And MEM[ j ] and MEM[ i No overlap, MEM[ j ] is the first j A MEM segment, MEM[ i ] is the first i MEM segments, len[ i ] indicates the first i The length of each MEM, cov[ i The `[]` represents the number of sequences it covers. If a conflict is found (i.e., two MEMs have overlapping positions on the same sequence), no state transition occurs. Finally, the backtracking of the `dp` array and the forward pointer `prev[` are used to determine the sequence number. i This method restores the globally optimal set of conflict-free MEMs. The time complexity of this method is O(...). k ²), suitable for medium-sized data.

[0043] The exact mode involves independently running dynamic programming for each maximum exact match, using the maximum fragment length as the objective function, eliminating overlaps sequence by sequence, and retaining the sequence that satisfies the minimum number of covering sequences. MThe maximum exact match, specifically, in exact mode, performs dynamic programming independently for each sequence, using the maximum segment length as the objective function, and eliminates overlaps sequence by sequence. For a given sequence, suppose there exists... m Given sorted MEMs, define the state dpp[ i [For the first] i The optimal solution ending with MEM has the following state transition:

[0044] dpp[ i ] = max(dpp[ j ] + len[ i ]); in, j < i And MEM[ j ].end <MEM[ i ].start.

[0045] That is, a combination of two MEMs is only considered if they have no positional conflict on the sequence. Finally, the optimal subset of MEMs on each sequence is selected, merged to form a global MEM set, and retains those that meet the preset minimum number of covered sequences. M MEMs. The time complexity of this method is O( d · k ²), where d This represents the total number of sequences and is suitable for comparison scenarios with many conflicts and where strict control of local overlap is required.

[0046] Both modes record paths using dynamic programming tables and predecessor pointers, supporting rapid backtracking and generating high-quality MEM sets that can be used for subsequent chain expansion and sequence segmentation. The fast mode focuses more on optimizing the overall area and is suitable for rapid preprocessing or high-density MEM regions; the precise mode is more suitable for low similarity or complex conflict scenarios, ensuring clean and non-redundant alignment partitions.

[0047] S103, expand each maximum exact match fragment in the maximum exact match fragment set to form a chain.

[0048] In one embodiment, each maximum exact match fragment in the maximum exact match fragment set is expanded to form a chain, including: expanding each maximum exact match fragment in the maximum exact match fragment set to all sequences using the SSW algorithm to form a chain.

[0049] Specifically, S3-1: Use the seed region of each MEM as a query sequence, and extract the corresponding positions before and after it. delta Length (e.g.) deltaThe fragments with a value of 20 are used as the expansion window; at the same time, the sequence fragments that are not covered by MEMs (i.e., the blank areas in the reference sequence) are extracted as the reference sequence.

[0050] S3-2: Utilizing the vectorized alignment core of the SSW algorithm (Single Instruction, Multiple Data, SIMD instruction acceleration), local alignment is performed between each MEM window and the reference sequence. A matching score of +2, a mismatch penalty of -3, a gap start penalty of -5, and a gap extension penalty of -1 are set. The optimal alignment path is calculated through dynamic programming, outputting a precise extension termination point. For example, when the MEM length is 8, if the alignment score to the left reaches a threshold (e.g., ≥10), the extension start point is moved to a new position; if the alignment to the right does not reach the threshold, it terminates at the end of the original MEM. All extended MEMs form a chain structure, recording the anchoring intervals of each sequence. The chain requires consistent spans (same starting offset) between sequences to maintain topological matching. The SSW algorithm's extension process is as follows: Figure 2 As shown.

[0051] S3-3: For sequence segments that cannot be expanded (such as alignment interruptions caused by repetitive regions), mark them as non-anchor intervals, temporarily store them in a buffer, and record their relative positions for processing in subsequent stages.

[0052] S104: The long sequence is segmented by the chain to obtain multiple subsequence sets; a third-party multi-sequence alignment software is called to align the multiple subsequence sets, and the alignment results of each subsequence set are spliced ​​together to obtain the multi-sequence alignment result.

[0053] In one embodiment, third-party multiple sequence alignment software is invoked to align multiple subsequence sets, including: using MAFFT's G-INS-i algorithm and HAlign3's fast mode to align multiple subsequence sets.

[0054] The sequence dataset is segmented using chains, dividing the sequences into multiple subsets. Third-party multi-sequence alignment software is then used to perform parallel alignment of these subsets to accelerate computation. Specifically, this includes the following sub-steps:

[0055] S4-1: Divide each sequence into multiple consecutive segments using the boundaries of the chain as dividing points. For example, if a sequence is divided into three regions [0-150], [151-300], and [301-500] by the chain, then these three regions constitute an independent set of subsequences.

[0056] S4-2: The G-INS-i algorithm of MAFFT (suitable for high similarity regions) and the fast mode of HAlign3 (suitable for low similarity regions) are adopted. Subsequence sets are submitted to the alignment engine in parallel through a multi-threaded pool, and a timeout threshold (e.g., 60 seconds / task) is set to prevent resource waste. The alignment results are saved as intermediate files.

[0057] S4-3: For misaligned subsets, trigger a recursive splitting strategy: if the subsequence length exceeds... If the length is less than 0.00bp, then dynamic window detection (window size w=50, step size s=10) is used again to find locally conservative fragments for secondary cutting, and a local alignment tool (such as MUSCLE) is called to reprocess it; if the length is insufficient... If the region is 0bp, it is marked as an unmatched region and proceeds to the subsequent sequence summary alignment stage.

[0058] In one embodiment, a profile alignment algorithm based on Fast Fourier Transform and K-Band constraints is used to align the maximum exact match segment that cannot be expanded. The alignment results of each subsequence set are concatenated with the alignment results of the maximum exact match segment that cannot be expanded to obtain a multi-sequence alignment result. The specific steps include the following: S5-1: Calculate the length of the unaligned region for each sequence and sort them in ascending order to prioritize the processing of shorter segments and reduce accumulated errors. For example, when the unmatched region length distribution is [15, 30, 5], prioritize processing segments of 5 bp.

[0059] S5-2: Profile alignment algorithm accelerated by Fast Fourier Transform (FFT): The aligned regions are constructed into a site-specific scoring matrix (PSSM), and the corresponding profile is generated (containing conservative weights for each position); for unmatched segments, FFT convolution is used to calculate their similarity score with the profile (time complexity O(n). d log d ), compared to traditional O( d ²) Significant optimization. Simultaneously, K-Band constraints (e.g., K=10) are introduced to limit the dynamic programming search range to comparisons only within the diagonal ±K stripes, further reducing computational cost (saving 80% of memory compared to ordinary algorithms), such as... Figure 3 As shown, Figure 3 This is a schematic diagram of the K-band process in dynamic programming. For example, for an unmatched region of length 50, a spectral signal is generated through FFT preprocessing and matched with the Profile spectrum to quickly identify highly similar segments.

[0060] S5-3: Combine all local alignment results into a complete global alignment: finally generate a complete multiple sequence alignment file in FASTA format.

[0061] In one embodiment, the present invention also provides a multi-sequence parallel alignment method based on maximum exact matching, such as... Figure 4 As shown, this embodiment includes three stages, specifically: Phase 1: Establish the suffix array and obtain MEMs.

[0062] Phase 2: Dynamic Programming Filtering MEMs.

[0063] Phase 3: Parallel multi-sequence alignment based on MEMs.

[0064] Parallel multi-sequence alignment includes SW extended MEMs, parallel alignment, Seq-Profile alignment, and fragment splicing.

[0065] In a specific embodiment, taking four DNA sequences as an example, the implementation process of the multi-sequence parallel alignment method based on maximum exact match provided by the present invention is illustrated: Input sequence: Sequence 1: ATCGTACGA Sequence 2: ATGGTACGCA Sequence 3: ATCGTAGGA Sequence 4: ATCGTAGTA Step S201: Sequence concatenation and suffix array construction.

[0066] The above sequences are concatenated into a long sequence, and then the gsaca-k algorithm is used to construct the suffix array (SA) and the longest common prefix (LCP) array. The concatenated long sequence is:

[0067] T = "ATCGTACGGA$ATGGTACGCA$ATCGTATGGA$ATCGTACGTA#" In this context, "$" is the separator. " is the end character.

[0068] Step S202: LCP-based MEM (maximum exact match) extraction.

[0069] Iterate through the LCP array, when LCP[ i When the number of occurrences is greater than or equal to 5 and originates from different sequences, extract candidate LCP intervals. For example:

[0070] LCP[4] = 6 corresponds to the prefix "ATCGTA", which appears in sequence 1, sequence 3 and sequence 4.

[0071] MEM1 = "ATCGTA", starting position [0, 0, 0], length 6.

[0072] Step S203: Use dynamic programming to filter non-overlapping ordered MEM sets.

[0073] With the goal of MEM "area = length × number of covered sequences", construct the DP table: MEM1(6x3=18), MEM2(5x2=10), MEM3(4x3=12).

[0074] Select the MEM set = {MEM1, MEM3}.

[0075] Step S204: Expand and form a chain based on the SSW algorithm Using MEM1 as the seed, perform SSW local alignment on each sequence with the following parameters: match: +2, mismatch: -3, initial penalty for gaps: -5, extension penalty: -1.

[0076] Comparison results: Sequence 1: ATCGTACGGA (score=28); Sequence 2: ATGGTACGCA (score=24); Sequence 3: ATCGTATGGA (score=25); Sequence 4: ATCGTACGTA (score=27).

[0077] Chain1 = [0–10].

[0078] Step S205: Segment the subsequence set and perform parallel alignment.

[0079] The long sequence is divided according to the boundary of Chain1 to obtain the subsequence sets: SubSet1 = [0–10], SubSet2 = [11–20], etc.

[0080] Call MAFFT G-INS-i or HAlign3 for parallel comparison, set the number of threads to 8, and the timeout threshold to 60s.

[0081] Step S206: Processing and splicing of misaligned segments.

[0082] For segments that have not been expanded, FFT+K-Band is used for fast comparison.

[0083] By stitching together all the partial results into a complete comparison, the final comparison result is as follows: Sequence 1: ATCGTACGA Sequence 2: AT-GGTACGCA Sequence 3: ATCGTAGGA Sequence 4: ATCGTAGTA The method provided by this invention primarily addresses the high complexity and low efficiency issues encountered when processing ultra-long DNA / RNA sequence data generated by third-generation sequencing, especially when the sequence length reaches or exceeds 1 Mbp (million base pairs). Traditional alignment methods often struggle to handle the massive data volume and complex sequences. Specifically, this invention, by combining a cutting strategy driven by multi-level suffix array indexing technology, enables precise localization and efficient processing of homologous regions across multiple sequences. The multi-level structure of the suffix array not only accelerates the sequence data indexing process but also improves matching accuracy. Redundant information and unnecessary calculations during the alignment process can be effectively removed, thereby improving alignment speed and accuracy. Furthermore, the multi-sequence alignment method of this invention enables efficient parallel processing, adapting to the needs of ultra-large-scale datasets, significantly reducing computational resource consumption, and shortening the time cost during the alignment process. This innovative parallel processing method can fully utilize modern high-performance computing resources, ensuring the efficiency and scalability of large-scale genomic data alignment. The optimized method of this invention has broad prospects in practical applications, particularly suitable for fields such as virus evolution tracking, species phylogenetic analysis, and genomics research. In the study of viral genome evolution, it can quickly and accurately identify genetic variations between viral strains, providing reliable data support for epidemic monitoring and early warning; in the analysis of species kinship, it can process and compare a large amount of genomic data from different species, helping scientists to more accurately infer the kinship and evolutionary history between species.

[0084] This invention improves the efficiency and accuracy of multiple sequence alignment by optimizing the chain finding and expansion process, and shows excellent performance, especially when dealing with large-scale, low-similarity and ultra-long sequence datasets.

[0085] The advantages of this invention include: 1. High efficiency: The suffix array and related structures are constructed using the gsaca-k algorithm, with a time complexity of O(N) and a space complexity of O(σ), which significantly reduces the computation and storage overhead of large-scale sequence data.

[0086] 2. Accuracy: By introducing the left expansion and dynamic programming filtering strategies of MEMs, the accuracy and orderliness of the maximum precise matching are ensured, conflicts and redundancy are avoided, and the reliability of the comparison results is improved.

[0087] 3. Scalability: The SSW algorithm is used to expand MEMs and form chains, which supports rapid expansion of local alignments, while retaining segments that cannot be expanded for subsequent processing, thus enhancing the adaptability of the algorithm.

[0088] 4. Parallel processing: By splitting the sequence dataset into chains and calling third-party multi-sequence alignment software (such as MAFFT and HAlign3) for parallel alignment, the computation process is significantly accelerated, making it suitable for ultra-long sequence datasets.

[0089] 5. Flexibility: It provides two dynamic programming strategies, fast mode and precise mode, which are optimized for global efficiency and local accuracy respectively to meet the needs of different scenarios.

[0090] 6. Optimize alignment quality: The FFT / K-Band-based Sequence-Profile alignment algorithm efficiently processes unmatched segments and generates complete multi-sequence alignment results through a splicing strategy, ensuring global consistency.

[0091] The inventive points of this invention include: 1. Suffix array construction based on gsaca-k algorithm: An improved induced sorting framework is adopted to efficiently generate suffix arrays, LCP arrays and document arrays, providing basic data structures for MEMs computation.

[0092] 2. Left expansion and dynamic programming filtering of MEMs: Left expansion ensures the non-scalability of MEMs, and dynamic programming strategy is combined to filter an ordered and conflict-free set of MEMs, thereby optimizing the accuracy and consistency of the comparison results.

[0093] 3. Chain extension of SSW algorithm: Using the vectorized alignment core of SSW algorithm, MEMs are extended to all sequences to form a chain structure, while retaining segments that cannot be extended for subsequent processing.

[0094] 4. Parallelized subsequence set alignment: By splitting the sequence dataset into chains, third-party multi-sequence alignment software is called for parallel processing, which significantly improves computational efficiency and is suitable for large-scale data.

[0095] 5. FFT / K-Band based Sequence-Profile alignment: The FFT-accelerated profile alignment algorithm, combined with K-Band constraints, efficiently processes unmatched segments, reduces computational complexity, and improves alignment quality.

[0096] 6. Recursive segmentation and splicing strategy: For subsequence sets that cannot be aligned, a recursive segmentation strategy is adopted for reprocessing, and all local alignment results are integrated through a splicing algorithm to generate a complete multi-sequence alignment file.

[0097] When applying the multi-sequence parallel alignment method based on maximum exact matching provided by this invention, it is not necessary to consider... Figure 1 The steps shown are executed in sequence. The specific execution order of each step can be determined as needed, and this invention does not impose any restrictions on it.

[0098] The above describes a multi-sequence parallel alignment method based on maximum exact matching, provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding multi-sequence parallel alignment device based on maximum exact matching, which includes: The module is used to concatenate strings from multiple sequences to obtain a long sequence, and uses a combination of bucketing and induced sorting to determine the suffix array and the longest common prefix array of the long sequence. The filtering module is used to find the maximum exact matching segment from the longest common prefix array based on the suffix array, and to filter the ordered and non-overlapping set of the maximum exact matching segments using dynamic programming techniques. The extension module is used to extend each maximum exact match fragment in the maximum exact match fragment set to all sequences, forming a chain; The alignment module is used to segment long sequences by chain to obtain multiple subsequence sets; it calls third-party multi-sequence alignment software to align multiple subsequence sets, and splices the alignment results of each subsequence set to obtain the multi-sequence alignment result.

[0099] Specific limitations regarding the maximum exact matching-based multi-sequence parallel alignment device can be found in the limitations of the maximum exact matching-based multi-sequence parallel alignment method described above, and will not be repeated here. Each module in the aforementioned maximum exact matching-based multi-sequence parallel alignment device can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the memory of a computer device in software form, so that the processor can call and execute the corresponding operations of each module.

[0100] The present invention also provides a computer-readable storage medium storing a computer program that can be used to execute the above-described... Figure 1 The provided method is a multi-sequence parallel alignment method based on maximum exact matching.

[0101] The present invention also provides Figure 5 The schematic diagram of the computer device shown is as follows: Figure 5 As shown, at the hardware level, this computer device includes a processor, internal bus, network interface, memory, and non-volatile memory, and may also include other hardware required for business operations. The processor reads the corresponding computer program from the non-volatile memory into memory and then executes it to achieve the above. Figure 1The provided method is a multi-sequence parallel alignment method based on maximum exact matching.

[0102] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the methods described above. Any references to memory, storage, databases, or other media used in the embodiments provided by this invention can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical storage, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0103] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this invention.

Claims

1. A multi-sequence parallel alignment method based on maximum exact matching, characterized in that, The method includes: Multiple sequences of strings are concatenated to obtain a long sequence. A combination of bucketing and induced sorting is used to determine the suffix array and the longest common prefix array of the long sequence. Based on the suffix array, find the maximum exact matching segment from the longest common prefix array, and use dynamic programming to filter the ordered and non-overlapping set of the maximum exact matching segments. Each maximum exact match fragment in the maximum exact match fragment set is expanded to all sequences to form a chain; Long sequences are split into multiple subsequence sets by using chains; A third-party multiple sequence alignment software is called to align multiple subsequence sets, and the alignment results of each subsequence set are spliced ​​together to obtain the multiple sequence alignment results.

2. The method according to claim 1, characterized in that, The long sequence contains separators "$" between different sequences and a terminator "" at the end of the long sequence. For other characters in the long sequence In lexicographical order, it is defined as .

3. The method according to claim 2, characterized in that, A combination of bucketing and induced sorting is used to determine the suffix array and the longest common prefix array of long sequences, including: Divide the long sequence into initial buckets based on the frequency of characters in the sequence; The algorithm recursively sorts the suffixes of special characters and refines the suffix array, longest common prefix array, and document array through two rounds of induced sorting.

4. The method according to claim 1, characterized in that, The suffix array includes the starting position of the suffix; based on the suffix array, find the maximum exact matching segment from the longest common prefix array, including: Iterate through the array of longest common prefixes to determine the longest common prefix interval; Get the initial position list, iterate through each longest common prefix interval, and add the starting position of the suffix of each longest common prefix interval to the position list; The left expansion checks whether the characters after subtracting the offset from the starting position of all suffixes are consistent. If they are consistent, the starting position of the suffix in the longest common prefix interval is expanded to the left by offset until it can no longer be expanded. The maximum exact matching segment is determined based on the expanded starting position and the length of the common prefix.

5. The method according to claim 4, characterized in that, The method further includes: If the length of the maximum exact match segment is greater than the preset sequence length, then the maximum exact match segment is discarded.

6. The method according to claim 1, characterized in that, Using dynamic programming techniques, select the ordered and non-overlapping set of the largest exact matching segments, including: For any maximum exact match segment, calculate the average position of each substring in the maximum exact match segment; If a duplicate substring appears in the same maximum exact match segment, delete the substring that deviates furthest from the average position; Sort all the most exact matching segments in ascending order of average position; By optimizing the conflict problem of all maximum exact matching segments using a dynamic programming strategy, an ordered and non-overlapping set of maximum exact matching segments is obtained.

7. The method according to claim 6, characterized in that, Dynamic programming strategies include fast mode and exact mode; fast mode involves using the maximum exact matching area as the objective function and filtering the set of conflict-free maximum exact matches that maximizes the total area through a single global dynamic programming operation. The exact mode involves independently running dynamic programming for each maximum exact match, using the maximum fragment length as the objective function, eliminating overlaps sequence by sequence, and retaining the sequence that satisfies the minimum number of covering sequences. M The maximum exact match.

8. The method according to claim 1, characterized in that, Each maximum exact match fragment in the maximum exact match fragment set is expanded to all sequences to form a chain, including: The SSW algorithm is used to extend each maximum exact match fragment in the maximum exact match fragment set to all sequences, forming a chain.

9. The method according to claim 1, characterized in that, Calling third-party multiple sequence alignment software to align multiple subsets of sequences, including: The G-INS-i algorithm of MAFFT and the fast mode of HAlign3 are used to compare multiple subsequence sets.

10. The method according to claim 1, characterized in that, The method further includes: The profile matching algorithm based on fast Fourier transform and K-band constraints is used to match the maximum exact matching segment that cannot be expanded. The alignment results of each subsequence set are concatenated with the alignment results of the largest exact match segment that cannot be expanded to obtain the multi-sequence alignment results.

Images