A low-memory multi-genome alignment and structural variation integration method for super large-scale closely related genome set

Abstract

This invention relates to a low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets. It solves the problems of existing multi-sequence / multi-genome alignment methods, which suffer from unacceptable time and space overhead on ultra-long sequences, large errors, and inefficiency in identification and integration. It includes S1, sequence input, output, and center sequence preprocessing; S2, direction determination and anchor point retrieval; S3, main strand construction, loop divide-and-conquer, and banded dynamic programming for fine alignment; and S4, structural difference strand identification, block output, and consistent integration. The advantages of this invention are: it can losslessly preserve and integrate structural variation fragments related to inversions and rearrangements during global alignment, ensuring that the output meets the "column consistency" requirements of downstream analysis while expressing strand direction and rearrangement information, and avoiding the extremely slow or even unworkable problems of traditional multiple merging processes in ultra-large file scenarios.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and comparative genomics, specifically to a low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome collections. Background Technology

[0002] With the continuous decline in sequencing costs, the number of genomic samples in public databases and clinical cohorts is growing rapidly. Unified alignment of thousands of human chromosome-level sequences, and even millions of viral whole genomes, is a crucial foundational step in population evolution, variant annotation, and precision medicine analysis. Existing multi-sequence / multi-genome alignment methods typically face three bottlenecks in ultra-large-scale scenarios: First, traditional full-matrix dynamic programming or complex graph models are prohibitively time- and space-consuming on ultra-long sequences; second, incremental merging relying on guide trees significantly amplifies errors and incurs high merging costs as the number of samples increases; third, structural variations and rearrangements such as inversions require the preservation of orientation and block boundary information in the results, but efficient identification and integration often contradict the goal of large-scale computation. Therefore, there is an urgent need for a multi-genome alignment technology that can scale to ultra-large-scale data, accommodate structural variation expression, and provide engineering-grade operability. Summary of the Invention

[0003] The purpose of this invention is to address the shortcomings of existing technologies for quickly constructing global multi-sequence (multi-genome) alignment results from closely related or same-species genome sets into closely related genome sets under conditions of significantly increased sample size and sequence length, while maintaining controllable memory usage. This invention proposes a low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets. This method can losslessly preserve and integrate structural variation fragments related to inversions and rearrangements during global alignment, ensuring that the output meets the downstream analysis requirements for "column consistency" while also expressing chain direction and rearrangement information. Simultaneously, it avoids the extremely slow or even unworkable problems of traditional multiplex merging processes in ultra-large file scenarios.

[0004] To achieve the above objectives, the present invention employs the following technical solution: a low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets, comprising the following steps:

[0005] S1. Sequence input, output, and center sequence preprocessing;

[0006] S2, Direction Judgment and Anchor Point Retrieval;

[0007] S3, main chain construction, loop divide and conquer and strip dynamic programming fine alignment;

[0008] S4, Structural Difference Chain Identification, Block Output and Consistency Integration.

[0009] In step S1, the input is a raw sequence file containing multiple nucleic acid sequences, and the output supports two formats:

[0010] ① One type is a commonly used multi-sequence alignment matrix file, which facilitates direct downstream statistics and visualization;

[0011] ② Multiple alignment result files organized in "blocks" express fragment orientation and rearrangement relationships, suitable for genome-wide alignment result management and exchange.

[0012] The center sequence preprocessing in step S1 specifically includes the following steps:

[0013] S11. For each original sequence Perform cleaning and indexing of records: characters that will interfere with subsequent alignments and coordinate mappings will be removed or normalized from the "working sequence" used for indexing and retrieval;

[0014] S12. Save the set of positions and counts of the characters to be removed in the original sequence as metadata;

[0015] S13, Central Sequence Under the default strategy, the sequence with the longest length after removing degenerate bases is selected;

[0016] In step S2, since homologous fragments in real genome data may appear in forward or reverse complementary directions, the orientation is determined before formal alignment. A rapid discrimination structure is constructed based on the short fragment characteristics of the center sequence, and the degree of sharing between the sequence to be aligned and the center sequence in both forward and reverse cases is evaluated, so as to select a more consistent orientation to enter the subsequent alignment process.

[0017] The method used is a fast discrimination based on k-mer set similarity:

[0018] ① Construct a forward k-mer set F and a reverse complementary k-mer set R from the central sequence;

[0019] ② The number of k-mer hits in F and R for the alignment sequence Q is calculated using the following formula:

[0020]

[0021] like If the condition is met, choose the positive direction; otherwise, choose the negative direction.

[0022] In the aforementioned low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome collections, during the anchor point retrieval stage, candidate homologous regions and their coordinate information are obtained by indexing the central sequence and quickly searching for common fragments between two sequences; or, mature genome alignment tools are used to first obtain candidate homologous regions, and then path filtering and merging are performed. The goal of the anchor point retrieval stage is to find the central sequence... Comparison with the sequence to be compared Quickly identify shared common segments between the two sequences and output their coordinates and length triples in both sequences. ,in For the center sequence coordinates, For the coordinates of the sequence to be compared, The segment length is specified; specifically, the following steps are included:

[0023] S21. Set a minimum anchor point length threshold. Filter out matches with short and unstable anchor lengths;

[0024] S22. The central sequence is established based on the FM-index of Burrows–WheelerTransform;

[0025] S23. Given a query fragment P, the FM-index supports... Solving for the set of occurrence positions ;

[0026] In step S3, the main chain construction includes the following steps:

[0027] S31. First, perform consistency filtering, using "offset consistency" as an intuitive criterion: define an offset for anchor point j. Within a local window, priority is given to preserving elements with the main offset. More consistent candidates;

[0028] S32. Subsequently, the candidate intervals are sorted by endpoints and non-overlapping filtering is performed, which is written as: for any adjacent incoming anchor points satisfy Under the constraint of [the condition], a greedy strategy is used to prioritize longer segments to improve skeleton stability;

[0029] S33. Finally, apply the "coordinate monotonically consistent" constraint to the remaining candidates to obtain the final main chain: Let the main chain be an index sequence. ,Require Under this constraint, the longest increasing subsequence idea is applied to construct a chain, and the objective function is to maximize the total anchor length of the chain: ;

[0030] In step S3, the cyclic divide-and-conquer method specifically includes the following steps:

[0031] S34. First, quickly split the data at a larger threshold to locate the homologous skeleton;

[0032] S35. Further refine within the sub-intervals until the interval length is within a controllable range, then switch to strip dynamic programming to complete fine alignment and gap filling.

[0033] That is: first set a larger minimum common segment threshold. Only retain Anchor points are used to quickly locate the main matching skeleton and establish a coarse-grained main chain. The main chain divides the two sequences into several corresponding sub-intervals, which can be formally represented as the r-th sub-interval. and Its boundary is induced by the endpoints of adjacent anchor points; then the threshold is reduced to within each subinterval. Then, retrieve denser anchor points and build finer sub-main chains. Iterate in this way until the interval length drops to the preset upper limit. Then switch to the fine alignment module.

[0034] In step S3, fine alignment employs strip dynamic programming, specifically including:

[0035] The banded regions are constructed around the aligned diagonal inferred from the main chain, with a length difference. Set baseline offset, only calculate if it meets the requirements. The state point; using affine gap penalty, let the cost of the gap of length t be... Where d is the penalty for opening a gap and e is the penalty for extending a gap;

[0036] The bandwidth is expanded incrementally or by doubling. Starting with a small expansion, it is gradually increased if the score or path indicates insufficient bandwidth, until the benefit no longer increases. A three-state recursion is used, with three matrices representing the optimal score ending in different "terminal states." The matching / mismatch scoring function is defined as follows: ,Right now This represents the score for aligning the i-th character of A with the j-th character of B;

[0037] definition:

[0038] M[i,j]: The best score when aligned to (i,j) and ending with "character-to-character"; I[i,j]: Aligned to (i,j) and ending with "inserting a gap in B";

[0039] D[i,j]: Aligned to (i,j) and ending with "Insert a notch in A";

[0040] Then, the strip dynamic programming performs the following recursion on the allowed in-band point (i,j):

[0041]

[0042] Where: I and D correspond to two different gap directions;

[0043] When transitioning from the gap state of M or the opponent to the current gap state, deduct the gap opening penalty d;

[0044] When extending continuously in the same gap direction, after deducting the extension penalty e, the final global alignment score is calculated in a strip-like manner, which relies on "recursion only for the state within the strip": for a given i, enumeration satisfies:

[0045] The value of j in the text.

[0046] In the aforementioned low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets, when the output is a traditional multiple sequence alignment matrix, star-shaped integration follows the consistency principle of "one-time insertion, permanent insertion": after each non-center sequence is aligned with the center sequence, a set of gap insertion patterns relative to the center sequence coordinates is generated; these gap patterns are then merged on the center sequence coordinate axis and synchronously propagated to all sequences, forming a multi-sequence alignment matrix with a unified column set. Since degenerate bases (e.g., N) may be removed during the preprocessing stage, the final output needs to be restored to the original coordinate system. The restoration method is to use the correspondence between "original sequence coordinates - alignment column coordinates" to insert the removed characters back to the vicinity of the correct alignment column, and coordinate with the gap propagation rules to avoid misalignment.

[0047] In the aforementioned low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets, step S4 specifically includes the following steps:

[0048] S41. In the block-organized output mode, the main chain is built in both the forward and reverse directions.

[0049] S42. Aggregate and prune necessary overlaps of inversion or rearrangement candidates formed by reverse consistency segments.

[0050] S43. Construct a candidate chain of structural differences by using candidate segments that were eliminated during the main chain screening process but still have local consistency characteristics.

[0051] S44. Joint screening is performed based on principles such as length, fragment spacing, coverage ratio, and consistency of local length between two sequences to retain more reliable structurally different chains.

[0052] S45. For the retained structural difference chains, use strip dynamic programming to fill the gaps between adjacent anchor points, thereby forming a writable structural difference alignment block.

[0053] Compared with the prior art, the advantages of the present invention are as follows:

[0054] I. This invention designs an overall workflow for ultra-large-scale sequence sets, which can maintain runnability even with a significant increase in sample size and sequence length. It also keeps memory usage within a controllable range through a memory-friendly implementation strategy, making it more suitable for actual computing scenarios involving millions of sequences or large-scale genomic data.

[0055] Second, this invention employs a combined strategy of "anchor point divide and conquer + strip dynamic programming": first, the sequence is segmented using anchor point alignment and the focus is on the regions of difference; then, precise alignment is completed within a limited bandwidth, thereby reducing the large amount of invalid table-filling computation required by traditional full-matrix dynamic programming and significantly improving the overall alignment efficiency. Strip dynamic programming is considered to be an effective way to reduce computational and spatial overhead in long sequence alignment.

[0056] Third, this invention introduces a direction prediction mechanism during the alignment process and combines it with the processing flow of forward and reverse main chains and inversion clusters, so that when there are structural variations such as inversions in the sequence, it can still maintain a stable anchor chain construction and alignment block expression ability, thereby improving the recognition and characterization effect of rearrangement events such as inversions.

[0057] Fourth, this invention uses a structural variation chain screening and lossless integration module to express the structural variation-related alignment blocks and global alignment results in a consistent manner. This avoids the performance bottleneck that may occur in the final merging stage and reduces the risk that structural variation information will be weakened or lost during multi-sequence merging, making the output results more conducive to subsequent structural variation analysis and annotation.

[0058] V. This invention simultaneously supports outputting multiple sequence alignment results in FASTA format and multiple genome alignment results in MAF format. MAF can express multiple alignments at the whole genome scale and supports information such as positive and negative strands, which is convenient for downstream evolutionary analysis, variant annotation and visualization. Attached Figure Description

[0059] Figure 1 This is a flowchart of the method of the present invention;

[0060] Figure 2 This is a comparison diagram of the effects of the present invention with those of other methods; Detailed Implementation

[0061] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0062] like Figure 1-2As shown, a low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets includes the following steps:

[0063] S1. Sequence input, output, and center sequence preprocessing;

[0064] S2, Direction Judgment and Anchor Point Retrieval;

[0065] S3, main chain construction, loop divide and conquer and strip dynamic programming fine alignment;

[0066] S4, Structural Difference Chain Identification, Block Output and Consistency Integration.

[0067] In step S1, the input is a raw sequence file containing multiple nucleic acid sequences, and the output supports two formats:

[0068] ① One type is a commonly used multi-sequence alignment matrix file, which facilitates direct downstream statistics and visualization;

[0069] ② Multiple alignment result files organized in "blocks" express fragment orientation and rearrangement relationships, and are suitable for the management and exchange of alignment results at the whole genome level;

[0070] This block format is commonly used in multiple sequence alignment processing and allows recording additional information such as the direction of the source chain.

[0071] Specifically, the center sequence preprocessing in step S1 includes the following steps:

[0072] S11. For each original sequence Perform cleaning and indexing of records: characters that will interfere with subsequent alignment and coordinate mapping (such as degenerate base N or other IUPAC degenerate symbols, gap symbols, etc.) will be removed or normalized from the "working sequence" used for indexing and retrieval;

[0073] S12. Save the set of positions and counts of the characters to be removed in the original sequence as metadata;

[0074] S13, Central Sequence Under the default strategy, the sequence with the longest length after removing degenerate bases is selected;

[0075] From an information content perspective, it can be regarded as the reference that best covers the homologous skeleton, thereby improving the stability of subsequent anchor point retrieval and chain filtering; at the same time, it allows users to explicitly specify the center sequence to meet specific reference preferences.

[0076] First, the input data is preprocessed to remove degenerate bases, gap symbols, and other characters that may affect subsequent alignment and coordinate mapping. At the same time, the position and number of degenerate bases in the original sequence are recorded as the basis for subsequent restoration of the original coordinate system.

[0077] By default, the center sequence is selected as the longest sequence after removing degenerate bases. After global alignment, the degenerate bases are inserted back according to the index mapping information saved in the preprocessing stage, based on the correspondence between the original sequence coordinates and the alignment coordinates, to ensure that the output coordinates are consistent with the original sequence. At the same time, through ordered writing and redundant column control, a large number of meaningless empty columns are avoided when degenerate bases appear in the same position in multiple sequences.

[0078] In step S2, since homologous fragments in real genome data may appear in either forward or reverse complementary directions, orientation determination is performed before formal alignment. This avoids costly anchor point retrieval and dynamic programming in the wrong direction, thereby improving the robustness and efficiency of large-scale alignment. A rapid discrimination structure is constructed based on the short fragment characteristics of the center sequence, and the degree of sharing between the sequence to be aligned and the center sequence in both forward and reverse orientations is evaluated. A more consistent orientation is then selected for subsequent alignment processes.

[0079] The method used is a fast discrimination based on k-mer set similarity:

[0080] ① Construct a forward k-mer set F and a reverse complementary k-mer set R from the central sequence;

[0081] ② The number of k-mer hits in F and R for the alignment sequence Q is calculated using the following formula:

[0082]

[0083] like If the condition is met, choose the positive direction; otherwise, choose the negative direction.

[0084] To control space consumption and accelerate queries on large-scale data, sets F and R are represented using approximate set structures such as Bloom filters, allowing for controllable false positives while significantly reducing memory usage and increasing throughput. Direction selection is based on the relative hit difference between the two directions, which is generally more robust to a small number of false positives. This type of direction inference mechanism based on k-mer and Bloom filters has been used in related implementations for targeted pre-screening of large-scale multi-genome alignments.

[0085] Furthermore, in the anchor point retrieval stage, candidate homologous regions and their coordinates are obtained by indexing the central sequence and quickly searching for common fragments between two sequences; or, mature genome alignment tools are used to first obtain candidate homologous regions, and then path filtering and merging are performed; thereby reducing the complexity of subsequent main chain construction and structural difference identification. Both of these methods follow the principle of "coarse localization first, then fine alignment." The goal of the anchor point retrieval stage is to find the central sequence... Comparison with the sequence to be compared Quickly identify shared common segments between the two sequences and output their coordinates and length triples in both sequences. ,in For the center sequence coordinates, For the coordinates of the sequence to be compared, The segment length is specified; specifically, the following steps are included:

[0086] S21. To limit the number of candidates and ensure retrieval speed, a minimum anchor length threshold is set. Filter out matches with short and unstable anchor lengths;

[0087] S22. In long sequence scenarios, in order to support fast substring queries, the central sequence is built with an FM-index based on Burrows–WheelerTransform;

[0088] S23, Given a query fragment P (e.g. from...) (The substring of length k is extracted from the middle), FM-index supports in Solving for the set of occurrence positions .

[0089] The retrieval complexity is usually linearly related to (|P|) (often written as O(|P|)), the actual constant is related to the sampling strategy), while the index storage is close to the compressed size.

[0090] The output of the anchor point retrieval phase can therefore be viewed as mapping from a large number of query fragments P) to several The collection provides input for the subsequent construction of the main chain.

[0091] For the retrieved common segments or candidate homologous intervals, firstly, consistency filtering is used to reduce duplicate matching and many-to-many interference; then, non-overlapping screening is performed on the candidate intervals, and further, a chain screening strategy with monotonically consistent coordinates is adopted to obtain the main chain anchor point combination that is consistent with the coordinate order of the center sequence.

[0092] In detail, anchor point sets often exhibit the phenomenon of "high coverage, small offset, but significant overlap," and one-to-many and many-to-many candidate matches are easily generated in repetitive sequence regions, leading to conflicts between different anchor points. Therefore, it is necessary to construct a "main chain" to represent the most consistent homologous skeleton of two sequences under the assumption of co-directional and colinearity. In step S3, the main chain construction includes the following steps:

[0093] S31. First, perform consistency filtering, using "offset consistency" as an intuitive criterion: define an offset for anchor point j. Within a local window, priority is given to preserving elements with the main offset. More consistent candidates; for example, retaining Smaller and Larger anchor points thus converge one-to-many candidates into an approximate one-to-one relationship, reducing discrete jumps introduced by repeating sequences;

[0094] S32. Subsequently, the candidate intervals are sorted by endpoints and non-overlapping filtering is performed, which is written as: for any adjacent incoming anchor points satisfy Under the constraint of [the condition], a greedy strategy is used to prioritize longer segments to improve skeleton stability;

[0095] S33. Finally, apply the "coordinate monotonically consistent" constraint to the remaining candidates to obtain the final main chain: Let the main chain be an index sequence. ,Require Under this constraint, the longest increasing subsequence idea is applied to construct a chain, and the objective function is to maximize the total anchor length of the chain: ;

[0096] When sequence lengths reach the level of chromosomes or the entire genome, performing full matrix dynamic programming directly over large differential intervals between anchor points results in time and space overhead that increases quadratically with the interval length. This can typically be roughly summarized as... The scale is too large to handle. Therefore, the method uses a cyclic divide-and-conquer approach to break down the large problem into many controllable smaller problems, and gradually refines the anchor point framework through threshold increments.

[0097] In step S3, the cyclic divide-and-conquer method specifically includes the following steps:

[0098] S34. First, quickly split the data at a larger threshold to locate the homologous skeleton;

[0099] S35. Further refine within the sub-intervals until the interval length is within a controllable range, then switch to strip dynamic programming to complete fine alignment and gap filling.

[0100] That is: first set a larger minimum common segment threshold. Only retain Anchor points are used to quickly locate the main matching skeleton and establish a coarse-grained main chain. The main chain divides the two sequences into several corresponding sub-intervals, which can be formally represented as the r-th sub-interval. and Its boundary is induced by the endpoints of adjacent anchor points; then the threshold is reduced to within each subinterval. Then, retrieve denser anchor points and build finer sub-main chains. Iterate in this way until the interval length drops to the preset upper limit. Then switch to the fine alignment module.

[0101] In step S3, fine alignment employs strip dynamic programming. The idea is that, given two highly similar pairs of nodes or nodes already constrained by the main chain to "near the diagonal," only the strip region of the DP matrix closest to the diagonal is calculated, thus reducing the complexity from... Down to , where n is the interval length and k is the bandwidth. Specifically, it includes:

[0102] The banded regions are constructed around the aligned diagonal inferred from the main chain, with a length difference. Set baseline offset, only calculate if it meets the requirements. The state point; to better reflect the biological insertion notch behavior, an affine notch penalty is adopted, where the cost of a notch of length t is... Where d is the penalty for opening a gap and e is the penalty for extending a gap;

[0103] This type of penalty scoring model is widely used in sequence alignment to avoid generating a large number of fragmented gaps. Meanwhile, a gradual or doubling strategy is employed for bandwidth, starting with a small expansion and progressively increasing it if the score or path indicates insufficient bandwidth, until the benefit no longer increases, thus achieving a trade-off between accuracy and overhead. To efficiently implement affine penalty scoring in dynamic programming, a three-state recursion is used, with three matrices representing the optimal score ending at different "terminal states." Let the matching / mismatch scoring function be... ,Right now This represents the score for aligning the i-th character of A with the j-th character of B;

[0104] definition:

[0105] M[i,j]: The best score when aligned to (i,j) and ending with "character-to-character"; I[i,j]: Aligned to (i,j) and ending with "inserting a gap in B";

[0106] D[i,j]: Aligned to (i,j) and ending with "Insert a notch in A";

[0107] Then, the strip dynamic programming performs the following recursion on the allowed in-band point (i,j):

[0108]

[0109] Where: I and D correspond to two different gap directions;

[0110] When transitioning from the gap state of M or the opponent to the current gap state, deduct the gap opening penalty d;

[0111] When extending continuously in the same gap direction, after deducting the extension penalty e, the key to calculating the final global alignment score in a strip is "recursion only for the state within the strip": for a given i, enumerate the following:

[0112] The value of j in the text.

[0113] Preferably, when the output is a traditional multiple sequence alignment matrix, star-shaped integration follows the consistency principle of "one-time insertion, permanent insertion": after each non-center sequence is aligned with the center sequence, a set of gap insertion patterns relative to the coordinates of the center sequence is generated; these gap patterns are merged on the coordinate axis of the center sequence and synchronously propagated to all sequences, forming a multiple sequence alignment matrix with a unified column set. Since degenerate bases (e.g., N) may be removed during the preprocessing stage, the final output needs to be restored to the original coordinate system. The restoration method is to use the correspondence between "original sequence coordinates - alignment column coordinates" to insert the removed characters back to the vicinity of the correct alignment column, and coordinate with the gap propagation rules to avoid misalignment.

[0114] The state outside the band is considered as (Unreachable), thus automatically preventing the optimal path from going out of bounds. Each row calculates at most approximately 2k+1 cells, with a total number of cells of approximately Complexity metric level .

[0115] Strip dynamic programming is a fundamental method in sequence alignment. By adopting a strip strategy that "limits the computational region", invalid table filling can be reduced, and time and space overhead can be lowered, making it more suitable for long sequences and large-scale tasks.

[0116] In terms of implementation, the system supports the use of common gap penalty models and allows parameter configuration. To further reduce memory usage and improve stability, engineering techniques such as "saving only necessary intermediate states", "compressing and storing backtracking information" and "uniformly pre-allocating memory" can be adopted to ensure stable operation under large-scale input.

[0117] In detail, step S4 includes the following steps:

[0118] S41. In the block-organized output mode, the main chain is built in both the forward and reverse directions.

[0119] S42. Aggregate and prune necessary overlaps of inversion or rearrangement candidates formed by reverse consistency segments in order to more accurately express the structural differences caused by changes in sequence orientation.

[0120] S43. Construct a candidate chain of structural differences by using candidate segments that were eliminated during the main chain screening process but still have local consistency characteristics.

[0121] S44. Joint screening is performed based on principles such as length, fragment spacing, coverage ratio, and consistency of local length between two sequences to retain more reliable structurally different chains.

[0122] S45. For the retained structural difference chains, use strip dynamic programming to fill the gaps between adjacent anchor points, thereby forming a writable structural difference alignment block.

[0123] To achieve a "consistent expression" between the structural difference blocks and the global alignment results, the blocks are sorted according to the center sequence coordinates, and the global alignment matrix is ​​segmented based on the block boundaries. The corresponding columns of the structural difference blocks are inserted or replaced, while the alignment columns that are duplicated are deleted to avoid duplicate expression and subsequent analysis deviations.

[0124] For scenarios involving merging multiple files or extremely large results, a hierarchical parallel merging and streaming read / write strategy is adopted to minimize the long-term resident memory of result files, thereby improving overall throughput and scalability.

[0125] The role of the co-linear main chain can be simply understood as follows: first, string together the most obvious and positionally consistent homologous segments from two sequences to form a global skeleton; this provides a stable reference for subsequent gap filling or merging of multi-sequence results. Highly consistent segments obtained from reverse retrieval, if exhibiting the characteristic of "opposite directions but good internal matching" in coordinate relationships, are likely to correspond to inversion candidates. The system will cluster and merge these candidate segments according to their mutual position and overlap, and prune overlapping edges when necessary to prevent the same inversion from being repeatedly expressed or having unclear boundaries. The intuitive purpose of this is to ensure that the output retains both the overall co-linear main chain and clearly presents structural differences such as "local direction reversals" as independent blocks.

[0126] In summary, the principle of this embodiment is as follows: global alignment is constructed using a "star-shaped organization" approach: a central sequence is selected as the anchor reference, so that each non-central sequence only needs to complete one double-sequence alignment with the central sequence, and then they are synthesized into a global multiple alignment according to a unified insertion rule. To ensure that this framework remains engineering-feasible for ultra-long and ultra-large-scale data, this invention combines several key steps:

[0127] First, a direction determination mechanism for closely related genomes is set up: a fast discrimination structure is constructed using a short fragment feature set of the central sequence. The degree of sharing with the central sequence is calculated for the sequence to be aligned in both forward and reverse orientations. Based on this, a more consistent orientation is selected for subsequent alignment, thereby handling inversion orientation differences without performing expensive global alignment. Such approximate set structures are often used in bioinformatics for efficient set discrimination and retrieval.

[0128] Secondly, anchor-driven homology interval retrieval and main chain construction are introduced: an index structure that supports fast common fragment retrieval is established for the central sequence to obtain candidate anchors, or in an alternative implementation, candidate homology intervals are generated first using mature genome alignment tools and then screened; the candidate anchors are sequentially subjected to consistency filtering and non-overlapping selection, and the idea of ​​"longest increasing subsequence" is used to eliminate crossover and disordered intervals to obtain a main chain anchor combination that is consistent with the coordinate order of the central sequence;

[0129] Furthermore, a cyclic divide-and-conquer strategy is adopted for ultra-long sequences: first, a skeleton is quickly established and the difference interval is shortened using stricter anchor conditions, and then the conditions are gradually relaxed within the difference interval to continue decomposition until the length of the sub-intervals enters a controllable range, after which fine-grained comparison is performed.

[0130] Finally, strip dynamic programming is used to fill and align the difference intervals, and path information compression and unified memory management are combined to achieve low memory backtracking. Strip dynamic programming takes advantage of the characteristic that "the optimal path of similar sequences is often close to the main diagonal", which can significantly reduce invalid computation and reduce time and space overhead.

[0131] In the multi-genome output mode, this invention further identifies structural variation-related chain fragments from the anchor set and forms structural variation alignment blocks. These blocks are then inserted or replaced into the global alignment matrix without loss, resulting in a final result that simultaneously contains sequential homologous regions and structural variation blocks. The output can generate commonly used multiple sequence alignment matrix files, as well as multiple alignment result files organized in blocks and expressing strand orientation and rearrangement information. The block format (such as MAF) serves as a universal exchange format for whole-genome multiple alignment and can record information such as source sequence, coordinates, and strand orientation.

[0132] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A low memory multi-genome alignment and structural variant integration method for very large sets of closely related genomes, characterized in that, Includes the following steps: S1. Sequence input, output, and center sequence preprocessing; S2, Direction Judgment and Anchor Point Retrieval; S3, main chain construction, loop divide and conquer and strip dynamic programming fine alignment; S4, Structural Difference Chain Identification, Block Output and Consistency Integration.

2. The method of claim 1, wherein, In step S1, the input is a raw sequence file containing multiple nucleic acid sequences, and the output supports two formats: ① One type is a commonly used multi-sequence alignment matrix file, which facilitates direct downstream statistics and visualization; ② Multiple alignment result files organized in "blocks" express fragment orientation and rearrangement relationships, suitable for genome-wide alignment result management and exchange.

3. The method of claim 2, wherein, The center sequence preprocessing in step S1 specifically includes the following steps: S11. For each original sequence Perform cleaning and indexing of records: characters that will interfere with subsequent comparisons and coordinate mappings will be removed or normalized from the "working sequence" used for indexing and retrieval; S12. Save the set of positions and counts of the characters to be removed in the original sequence as metadata; S13, Central Sequence The sequence selected under the default strategy is the one with the longest length after removing degenerate bases.

4. The method for low-memory multi-genome alignment and structural variation integration for ultra-large-scale closely related genome sets according to claim 3, characterized in that, In step S2, since homologous fragments in real genome data may appear in forward or reverse complementary directions, the direction is determined before formal alignment; a rapid discrimination structure is constructed based on the short fragment characteristics of the center sequence, and the degree of sharing between the sequence to be aligned and the center sequence in both forward and reverse cases is evaluated, so as to select a more consistent direction to enter the subsequent alignment process. The method used is a fast discrimination based on k-mer set similarity: ① Construct a forward k-mer set F and a reverse complementary k-mer set R from the central sequence; ② The number of k-mer hits in F and R for the alignment sequence Q is calculated using the following formula: like If the condition is met, choose the positive direction; otherwise, choose the negative direction.

5. The method for low-memory multi-genome alignment and structural variation integration for ultra-large-scale closely related genome sets according to claim 4, characterized in that, In the anchor point retrieval stage, candidate homologous regions and their coordinates are obtained by indexing the central sequence and quickly searching for common fragments between two sequences; or, mature genome alignment tools are used to first obtain candidate homologous regions, and then path filtering and merging are performed; the goal of the anchor point retrieval stage is to find the central sequence Comparison with the sequence to be compared Quickly identify shared common segments between the two sequences and output their coordinates and length triples in both sequences. ,in For the center sequence coordinates, For the coordinates of the sequence to be compared, The segment length is specified; specifically, the following steps are included: S21. Set a minimum anchor point length threshold. Filter out matches with short and unstable anchor lengths; S22. The central sequence is established based on the FM-index of Burrows–WheelerTransform; S23. Given a query fragment P, the FM-index supports... Solving for the set of occurrence positions .

6. The method for low-memory multi-genome alignment and structural variation integration for ultra-large-scale closely related genome sets according to claim 5, characterized in that, In step S3, the main chain construction includes the following steps: S31. First, perform consistency filtering, using "offset consistency" as an intuitive criterion: define an offset for anchor point j. Within a local window, priority is given to preserving elements with the main offset. More consistent candidates; S32. Subsequently, the candidate intervals are sorted by endpoints and non-overlapping filtering is performed, which is written as: for any adjacent incoming anchor points satisfy Under the constraint of [the condition], a greedy strategy is used to prioritize longer segments to improve skeleton stability; S33. Finally, apply the "coordinate monotonically consistent" constraint to the remaining candidates to obtain the final main chain: Let the main chain be an index sequence. ,Require Under this constraint, the longest increasing subsequence idea is applied to construct a chain, and the objective function is to maximize the total anchor length of the chain: .

7. A low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets according to claim 6, characterized in that, In step S3, the cyclic divide-and-conquer method specifically includes the following steps: S34. First, quickly split the data at a larger threshold to locate the homologous skeleton; S35. Further refine within the sub-intervals until the interval length is within a controllable range, then switch to strip dynamic programming to complete fine alignment and gap filling. That is: first set a larger minimum common segment threshold. Only retain Anchor points are used to quickly locate the main matching skeleton and establish a coarse-grained main chain. The main chain divides the two sequences into several corresponding sub-intervals, which can be formally represented as the r-th sub-interval. and Its boundary is induced by the endpoints of adjacent anchor points; then the threshold is reduced to within each subinterval. Then, more dense anchor points are retrieved and finer sub-main chains are built. Iterate in this way until the interval length drops to the preset upper limit. Then switch to the fine alignment module.

8. The method for low-memory multi-genome alignment and structural variation integration for ultra-large-scale closely related genome sets according to claim 7, characterized in that, In step S3, the precise alignment using strip dynamic programming specifically includes: The strip-shaped regions are constructed around the aligned diagonal inferred from the main chain, with a length difference. Set baseline offset, only calculate if it meets the requirements. The state point; using affine gap penalty, let the cost of a gap of length t be... Where d is the opening gap penalty and e is the extension gap penalty; The bandwidth is expanded incrementally or by doubling. Starting with a small expansion, it is gradually increased if the score or path indicates insufficient bandwidth, until the benefit no longer increases. A three-state recursion is used, with three matrices representing the optimal score ending in different "terminal states." The matching / mismatch scoring function is defined as follows: ,Right now This represents the score for aligning the i-th character of A with the j-th character of B; definition: M[i,j]: The best score when aligned to (i,j) and ending with "character-to-character"; I[i,j]: The score when aligned to (i,j) and ending with "inserting a gap in B"; D[i,j]: Aligned to (i,j) and ending with "Insert a notch in A"; Then, the strip dynamic programming performs the following recursion on the allowed in-band point (i,j): Where: I and D correspond to two different gap directions; When transitioning from the gap state of M or the opponent to the current gap state, deduct the gap opening penalty d; When extending continuously in the same gap direction, after deducting the extension penalty e, the final global alignment score is calculated in a strip-like manner, which relies on "recursion only for the state within the strip": for a given i, enumeration satisfies: The value of j in the text.

9. A low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets according to claim 8, characterized in that, When the output is a traditional multiple sequence alignment matrix, star-structure integration follows the consistency principle of "one-time insertion, permanent insertion": after each non-center sequence is aligned with the center sequence, a set of gap insertion patterns relative to the coordinates of the center sequence is generated; these gap patterns are merged on the coordinate axis of the center sequence and synchronously propagated to all sequences to form a multiple sequence alignment matrix with a unified column set; since degenerate bases are removed in the preprocessing stage, the final output needs to be restored according to the original coordinate system. The restoration method is to use the correspondence between "original sequence coordinates - alignment column coordinates" to insert the removed characters back to the vicinity of the correct alignment column, and coordinate with the gap propagation rules to avoid misalignment.

10. A low-memory multi-genome alignment and structural variation integration method for ultra-large-scale closely related genome sets according to claim 9, characterized in that, Step S4 specifically includes the following steps: S41. In the block-organized output mode, the main chain is built in both the forward and reverse directions. S42. Aggregate and prune necessary overlaps of inversion or rearrangement candidates formed by reverse consistency segments. S43. Construct a candidate chain of structural differences by using candidate segments that were eliminated during the main chain screening process but still have local consistency characteristics. S44. Joint screening is performed based on principles such as length, fragment spacing, coverage ratio, and consistency of local length between two sequences to retain more reliable structurally different chains. S45. For the retained structural difference chains, use strip dynamic programming to fill the gaps between adjacent anchor points, thereby forming a writable structural difference alignment block.

Images