A multi-genome alignment method that can be used for large-scale genomes

Abstract

The scheme discloses a multi-genome alignment method which can be used for large-scale genomes, belongs to the technical field of computer biology, and proposes a BWT-FM-LIS cyclic division and governance mode to search for common substrings among whole genomes, searches for as many common substrings and key paths among genomes as possible, shortens the length of difference substring which needs dynamic programming alignment, and further reduces the complexity of alignment, so that the scheme can be used for alignment of large-scale large genomes, is suitable for large-scale long sequence or ultralong genome data, and can identify various structural variations in genome alignment.

Description

Technical Field

[0001] This approach belongs to the field of computer biology technology and proposes a multi-genome alignment method that can be used for large-scale genomes. Background Technology

[0002] Currently, the most mature genome sequence alignment tools are double sequence alignment, which has been widely used in genomics research. Most tools employ a "seed-expansion" method, such as BLAST. First, short, gap-free sequence matches (seeds) are identified. Then, an improved Smith-Waterman algorithm is used to expand the matches from both ends of the seed sequences until the alignment score falls below a specified threshold. Seed types can be classified according to whether precision is required, whether continuous matching is required, or whether the match length is fixed. BLAT and STELLAR use k-mers as their precise matching seeds, LASTZ (an updated version of BLASTZ) uses spaced seeds, and LAST utilizes suffix arrays to find adaptive seeds of varying lengths. MUMmer and CHAOS use suffix trees and thread trees, respectively, to quickly find all precise continuous matching seeds with a certain minimum length. Adaptive, imprecise, or spaced seeds often have higher sensitivity. After expansion, adjacent, sequential, and oriented local alignments can be linked together by linking programs such as AXTCHAIN ​​to form larger alignments as output.

[0003] The development of multiple sequence alignment tools for genomes started later than that of double sequence alignment, and their enormous temporal and spatial complexity has limited their application. Currently developed genome multiple sequence alignment tools mainly include hierarchical and local alignment methods. Representative tools for hierarchical alignment include: M-GCAT, progressiveMauve, Mugsy, progressiveCactus, Parsnp, and FAME; Parsnp, MAVID, MLAGAN, SeqAn::T-Coffee, and PECAN, which can perform progressive global multiple sequence alignment; DIALIGN, which uses greedy global multiple sequence alignment; and FSA, which is based on a probabilistic statistical hidden Markov model.

[0004] Currently, the most powerful interspecific genome multiple sequence alignment software is ProgressiveCactus (hereinafter referred to as Cactus), which can align approximately 600 amniotic genomes. However, its limitation is that it requires a scientific guide tree as prior knowledge input and cannot handle characters compatible with Windows. The most powerful intra-population genome multiple sequence alignment software is Parsnp, which can align thousands of microbial genomes. However, its drawback is the "weakest link" effect, meaning it can only find regions common to all sequences, leading to a significant decrease in alignment quality as the number of sequences increases. Furthermore, it cannot handle data with low similarity or large length differences. Although Parsnp and Cactus have been successfully applied to tens of thousands of mitochondrial genomes, thousands of highly similar microbial genomes, and approximately 600 amniotic genomes, respectively, there is currently no software capable of aligning thousands of mammalian genomes or even larger-scale multi-genome alignments. Summary of the Invention

[0005] The purpose of this solution is to address the aforementioned issues by providing a multi-genome alignment method and system applicable to large-scale genomes. This method is suitable for large-scale long-sequence or ultra-long genome data and can identify various structural variations in genome alignment.

[0006] A multi-genome alignment method applicable to large-scale genomes, comprising:

[0007] S1. Perform data preprocessing on the genome sequences to be aligned and determine the center sequence;

[0008] S2. Establish a BWT index for the center sequence to obtain the BWT data structure, and build a prediction dictionary;

[0009] S3. Use a prediction dictionary to predict the direction of each non-center sequence;

[0010] S4. Based on the common substring length threshold, perform a forward common substring retrieval for each non-center sequence within the BWT data structure of the center sequence. Here, "forward" refers to the direction relative to the center sequence; the forward direction of the center sequence is considered the forward direction, and the reverse direction is considered the reverse direction. For example, if the direction of a non-center sequence is AB and the direction of the center sequence is BA, then a one-way BA retrieval is performed on the non-center sequence.

[0011] S5. Select the longest common substring combination as the main chain of the corresponding non-center sequence;

[0012] S6. Filter the main chain to obtain pairs of differing substrings. Determine if the length of the differing substrings satisfies dynamic programming. If yes, continue to S7. Otherwise, remove the main chain and lower the threshold for the length of the common substring, and repeat S4-S6. The intervals corresponding to the main chains of the non-center sequence and the center sequence are exactly the same and are aligned. There is no need to continue comparing them. Just proceed to the final result. Except for the first loop, in each subsequent execution of step S4, the retrieved center and non-center sequences are the sequences after removing the latest main chain.

[0013] S7. Use dynamic programming to perform a fine comparison of the differences in substring pairs;

[0014] S8. Output the results of multiple sequence alignment.

[0015] This scheme proposes a BWT-FM-LIS cyclic divide-and-conquer approach to retrieve common substrings across the entire genome, aiming to find as many common substrings and critical paths as possible, thereby shortening the length of the differential substrings that require dynamic planning for alignment. This reduces the complexity of the alignment, making the scheme applicable to large-scale genome alignments.

[0016] In the aforementioned multi-genome alignment method applicable to large-scale genomes, in step S1,

[0017] Data preprocessing includes: recording the positions of degenerate bases in binary form and removing them along with the gaps;

[0018] The center sequence can be determined by either selecting the longest sequence after removing degenerate bases, or by specifying it by the user.

[0019] In step S2, two hash dictionaries are built for the forward and reverse chains of the center sequence to obtain the prediction dictionary. The prediction dictionary includes the forward chain Kmer dictionary and the reverse chain Kmer dictionary of the center sequence.

[0020] Step S2 also includes storing the central sequence using a binary compression storage method;

[0021] In step S3, the number of Kmers shared by each non-central sequence is found in the forward chain Kmer dictionary and the reverse chain Kmer dictionary of the central sequence, and the direction of the non-central sequence is determined accordingly.

[0022] In the aforementioned multi-genome alignment method applicable to large-scale genomes, the forward common substring retrieval method in step S4 is as follows:

[0023] S1. Starting from the beginning position B_index of the non-center sequence, the initial len is 0. Continuously increment len ​​by 1 to search for all common substrings in the center sequence until there are no common substrings. Backtrack one step to get x common substrings under the current step length len.

[0024] S2. Determine whether len is greater than or equal to the common substring length threshold. If so, record the common substring information found under that len, and set B_index = B_index + len, then return to step S1 to repeat the search.

[0025] If not, then set B_index = B_index + 1 and return to step S1 to repeat the search;

[0026] S3. The above process continues to retrieve the sequence until the end of the non-center sequence.

[0027] ACGTAAACGT

[0028] ACGTGACGA

[0029] Taking the above sequence pair as an example, the top one is the center sequence and the bottom one is the non-center sequence. In this example, starting from A in the non-center sequence, the initial step length len is 0, and there are 4 common substrings. Then add 1 to make it 1. Find AC in the center sequence, and there are 2 common substrings. Then add 1 to make it 2. Find ACG, and there are 2 common substrings. Then add 1 to make it 3. Find ACG, and there are 2 common substrings. Then add 1 to make it 4. Find ACGT, and there are 2 common substrings. Then add 1 to make it 5. Find ACGT, and no common substrings can be found.

[0030] Then we can get 4 as the two common substrings.

[0031] Assuming the threshold is 3, then len = 4 is greater than the threshold 3, so the information of the two common substrings with len = 4 is recorded. Then, starting from the starting position = B_index + 4, with len initially set to 0, the above steps are continued, that is, starting from G in the non-center sequence with an initial step size of 0, and no common substring can be found.

[0032] ACGT

[0033] CCGT

[0034] In the example above, after C moves one step, it becomes CC. The center sequence has no common substring. Assuming the threshold is greater than 1, then len is less than the threshold. Starting from the second C, when len = 0, a common substring is found; when len = 1, a common substring CG is found; when len = 2, a common substring CGT is found; when len = 3, there is no common substring, so when len = 2, there is a common substring CGT.

[0035] In the aforementioned multi-genome alignment method applicable to large-scale genomes, the recorded common substring information is in a one-to-one or one-to-many pattern depending on the retrieval results.

[0036] The methods for finding the main chain in step S5 include:

[0037] S51. First, all common substrings retrieved for the corresponding non-center sequence are filtered into a one-to-one mode according to the principle of closest relative position;

[0038] S52. Find overlapping and intersecting common substrings. For each pair of overlapping and intersecting common substrings, retain the longest corresponding homologous interval in the non-center sequence and discard the other corresponding homologous intervals in the non-center sequence to filter out the longest homologous interval combination. The longest non-overlapping and non-intersecting common substring combination is the main chain.

[0039] Some common substrings are one-to-one, and some are one-to-many. Here, we uniformly select a one-to-one pattern, meaning a segment on a non-central sequence corresponds to a segment on a central sequence. The main strand is composed of common substrings, which are small, identical segments between two sequences. Therefore, for each non-central sequence's main strand, there is a corresponding one in the central sequence. We select the one-to-one pattern according to the principle of relative position. This ensures that the corresponding main strand is selected from the central sequence while minimizing the offset between the two sequences, resulting in more consistent subsequent alignments. Furthermore, this scheme utilizes the remaining common substrings of the main strand to find structural variation chains. We use common substrings from the entire genome to identify regular structural variation chains, and then fill in the gaps in these chains to obtain structural variation regions for alignment, thus realizing a scientific method for extracting structural variations.

[0040] In the multi-genome alignment method described above that can be used for large-scale genomes, step S7 employs K-band dynamic programming to perform fine-grained alignment of differential substring pairs:

[0041] First, initialize the state table;

[0042] Then, the score of the current row is calculated using affine penalty, and the path transition direction is filled into the state table based on the relationship between the scores of the current row and the previous row.

[0043] Finally, the optimal path is derived by backtracking from the state table.

[0044] Among the aforementioned multi-genome alignment methods applicable to large-scale genomes, Kband dynamic programming,

[0045] Affine penalty is used to differentiate between first and second insertions, so that the penalty for the first insertion is greater than that for the second insertion;

[0046] Save time and space costs by limiting the scope of dynamic programming: restrict dynamic programming to the vicinity of the diagonal, and only calculate the region near the diagonal; the baseline of the strip is determined by the difference in length between the two sequences, and based on the strip, expand to the upper right and lower left, starting from 1 and increasing in powers of 2, until the score of the next change is no higher than the score of the previous change.

[0047] Simplify the data structure by storing path information: During the calculation process, a separate information table is set up to store the path information. When backtracking, the best path is found directly by looking up the path information table. During the calculation process, the current three states are stored, and then each row is calculated to overwrite the previous row. The path transfer direction of each state is recorded throughout the process.

[0048] Instead of repeatedly allocating and freeing memory, allocate global space: allocate the maximum memory that may be used in a single allocation, and replace repeated allocation and freeing of memory with repeated read and write operations.

[0049] In the above-mentioned multi-genome alignment method applicable to large-scale genomes, in step S7, the differential substring pairs are finely aligned by star alignment to obtain the insertion information of each sequence.

[0050] Furthermore, in step S7 of this method, post-processing is performed on the results after fine comparison:

[0051] The degenerate bases removed during preprocessing are integrated into the final alignment;

[0052] The degenerate base insertion process is as follows: first, the length after insertion is calculated; then, a vector of the corresponding size is allocated; and the original string information and the insertion information are written into the vector in sequence to realize the insertion operation.

[0053] In the above-mentioned multi-genome alignment method applicable to large-scale genomes, when step S8 is required to output the structural variation alignment results (the user can choose to output multiple sequence alignment results (FSATA mode) or obtain two results simultaneously (MAF mode) as needed), in step S4, a bidirectional common substring search is performed on each non-center sequence to find possible inverted structural variations.

[0054] Step S5 will obtain a forward main chain composed of the longest common substrings in the forward direction and a reverse main chain composed of the longest common substrings in the reverse direction. The forward main chain and the reverse main chain will be merged to obtain the final main chain of the corresponding non-center sequence. The merging process includes:

[0055] Find the forward main chain and the reverse main chain; map the reverse main chain to the forward main chain; merge the forward main chain and the mapped reverse main chain.

[0056] Among the aforementioned multi-genome alignment methods applicable to large-scale genomes, the forward retrieval method in bidirectional searching is as follows:

[0057] Starting from the beginning of the non-center sequence, after a step length len, if x common substrings are found in the center sequence, then determine whether the length of len is greater than the common substring length threshold.

[0058] If so, the search results are recorded as common substring information, and then the above search is repeated from the position of step length + len.

[0059] If not, then continue the above search from the position of step size + 1;

[0060] The above process continues until the end of the non-center sequence;

[0061] The reverse retrieval method is:

[0062] Starting from the beginning of the non-center sequence, after a step length len, if x common substrings are found in the center sequence, then determine whether the length of len is greater than the common substring length threshold.

[0063] If so, the search results are recorded as common substring information, and then the above search is repeated from the position of step length + len.

[0064] If not, then continue the above search from the position of step size + 100;

[0065] The above process continues until the end of the non-center sequence;

[0066] In the aforementioned multi-genome alignment methods applicable to large-scale genomes, this method further includes:

[0067] The common substrings filtered out by the forward main chain and the reverse main chain are used as the data source for finding other structural variant chains;

[0068] Two search operations were performed on the data source to find the positive structural variation chain and the structural variation chain combined with the inversion;

[0069] Selective retention processing is applied to overlapping structural variant chains;

[0070] In step S7, the dynamic programming method is used to compare and fill the regions between common substrings of the structural mutation chain to obtain the entire aligned structural mutation alignment, and write it to the maf file;

[0071] Multiple sequence alignment results are written to a FASTA file;

[0072] MAF files are used to store structural variation alignment results, while FASTA files are used to store multiple sequence alignment results.

[0073] In step S8, the maf file and fasta file are merged to output the final multi-genome alignment results:

[0074] First, sort the blocks in the MAF file according to the center sequence index and then filter the sequences in the block whose alignment quality exceeds that of the original FASTA file.

[0075] Then, based on the start and end points of each sequence in the block as cutting points, the FASTA alignment matrix is ​​column-cut.

[0076] Finally, based on the star comparison, the blocks of the MAF file and the corresponding blocks after the FASTA file are cut are integrated, and the structural variations in the original FASTA file are deleted.

[0077] Before merging, the MAF files are sorted and the sorting results are saved. Then, each block is read in and merged, and the merged result is written out. The merging process adopts a layer-by-layer merging method to ensure that the merging within each level is independent of each other.

[0078] The above method enables the algorithm provided in this scheme to output structural variant MAF files, and finally merges the structural variant MAF files and the multiple sequence alignment result FASTA files to output the final MAF format multi-genome alignment results that integrate the structural variants.

[0079] The advantages of this solution are:

[0080] (1) This scheme innovates a method for retrieving common substrings between genomes using a cyclic divide-and-conquer approach. By combining algorithms such as BWT with linear time and space complexity, data compression, FM-index, and longest increasing subsequence (LIS) in the BWT-FM-LIS cyclic divide-and-conquer approach, it retrieves as many common substrings and critical paths between genomes as possible, shortens the length of the differential substrings that need to be dynamically planned for alignment, and thus reduces the complexity of the alignment so that this scheme can be used for large-scale genome alignment.

[0081] (2) This scheme introduces affine penalty, path storage, Kband constraint and other methods in large-scale genome multiple sequence alignment for dynamic planning of optimizing highly similar genomes. While improving the biological significance of the alignment results, it reduces the computation time and storage space in the dynamic planning process. In large-scale human whole genome multiple sequence alignment, the dynamic planning of the differential substrings accounts for the highest proportion of time and space complexity. This improvement can greatly reduce the computation time and storage space of the alignment.

[0082] (3) This scheme uses common substrings of the whole genome to identify regular structural variation chains, and then fills the chain holes to obtain structural variation region alignment, realizing a scientific method for extracting structural variations.

[0083] (4) This scheme proposes a method for integrating structural variations with MSA results, and adopts an innovative star alignment strategy that is more efficient than multiz to integrate structural variations into multiple sequence alignment results;

[0084] (5) The software based on this scheme is the first software capable of aligning thousands of human chromosome 1 genomes, greatly increasing the upper limit of the data scale that can be processed for multi-genome alignment, and enabling large-scale genome alignment tasks that no current software can accomplish. Moreover, this scheme also demonstrates sufficiently excellent quality for small-scale genome alignment tasks that other software can complete, while significantly reducing the consumption of time and space costs. Attached Figure Description

[0085] Figure 1 This is an overall flowchart of the implementation of the present invention;

[0086] Figure 2 This is a summary of the specific algorithm flow in the process of this invention, wherein,

[0087] A represents an example of compressing a DNA sequence to 1 / 4 of its original space;

[0088] B represents an example process of finding the main chain in the common substring;

[0089] C represents the process of merging the positive and negative main chains;

[0090] D is an example of searching for strictly structural variant chains;

[0091] E is a diagram of the affine penalty;

[0092] F is a specific example of K-band dynamic programming;

[0093] G represents compressed storage of the backtracking path table;

[0094] H represents a specific example of deleting N and inserting N into a 3-sequence array;

[0095] I indicates the specific process of integrating the MAF result file of structural variation into the FASTA result file of multiple sequence alignment;

[0096] Figure 3 The results show the comparative analysis of Cactus, Parsnp, and Halign-G on 24 sets of chromosome datasets from 18 individuals.

[0097] A represents the distribution of the M-scores of the three software programs;

[0098] B represents the time and memory consumption of the three software programs during the experiment;

[0099] Figure 4The analysis shown compares seven software programs—Cactus, Parsnp, Mugsy, Kalign, MLAGAN, PMauve, and Halign-3—with Halign-G's FSATA and MAF modes on five datasets.

[0100] A represents the distribution of the M-scores for the nine software programs;

[0101] B represents the time and memory consumption of the nine software programs during the experiment;

[0102] Figure 5 The image shows a comparison of the structural variations identified by Cactus, Parsnp, Mugsy, Pmauve, and Halign-G across nine simulated datasets. Detailed Implementation

[0103] The present solution will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0104] This embodiment presents a multi-genome alignment method applicable to large-scale genomes. The main features include a cyclic divide-and-conquer approach for retrieving common substrings across the entire genome, a dynamic programming method for optimizing highly similar genomes, a scientific method for extracting structural variations, and a method for integrating structural variations with MSA results. This embodiment uses the human chromosome 1 genome as an example to illustrate this approach.

[0105] See the overall flowchart Figure 1 The following are the specific steps:

[0106] S1. Data Preprocessing and Determining the Central Sequence

[0107] This solution uses the FASTA file format as the input data format. The input data can contain various degenerate bases and gap characters. The input format can be either a single FASTA file containing multiple sequences, or a specified folder, which the software will process through all files within that folder. This solution uses binary pairs.<index,number> The positions of degenerate bases are recorded in a format that is then removed along with any gaps. The sequence type (RNA or DNA) is determined based on its characteristics. Furthermore, this scheme converts all remaining nucleic acid characters to uppercase and writes them to an intermediate file for later retrieval from external storage rather than resident in memory, facilitating the processing of extremely large datasets. The longest sequence after removing degenerate bases contains the most valid information and is therefore selected as the center sequence. Of course, in practical applications, the center sequence can also be specified by the user.

[0108] Figure 1As shown, in the original sequence, the longest DNA sequence is the identified central sequence. Among the other sequences, the top few are forward non-central sequences, and the bottom one is a reverse non-central sequence. The black pentagram represents degenerate bases. The process surrounded by the dark arrow in the dashed box is the BWT-FM-LIS cycle divider.

[0109] S2. Build a BWT index and prediction dictionary for the center sequence.

[0110] Create a BWT index for the central sequence

[0111] The BWT index construction process employs the best available suffix array constructor—the parallelizable LibDivSufSort—which can construct a suffix array in linear time complexity. The constructor is a function, and the BWT index can be viewed as a suffix array; building the BWT index involves obtaining this suffix array through the constructor's function. The BWT data structure of the central sequence is obtained by building the BWT index for the central sequence.

[0112] Build a prediction dictionary for the center sequence

[0113] In multi-genome alignment, genome sequences from different species or individuals are typically processed, and these sequences may have different orientations. This scheme uses MurmurHash3 technology to build two hash dictionaries for the Kmer of the central sequence's forward and reverse chains—a forward chain Kmer dictionary and a reverse chain Kmer dictionary. Subsequently, the hash key will be used to perform a fast search in the two dictionaries.

[0114] Central sequence compressed storage

[0115] To conserve memory during alignment, this scheme employs binary compressed storage for the center sequence. This compressed storage method uses each byte to store 4 bases, with each base occupying 2 bits. By compactly storing multiple bases within a single byte (see...),... Figure 2 In the middle (A), memory usage can be significantly reduced (to one-quarter of the original), increasing the scale of comparable data.

[0116] S3. Search for common substrings

[0117] kmer predicts direction

[0118] Find the number of Kmers shared by each non-central sequence in the Kmer dictionaries of the forward and reverse chains of the central sequence. Determine whether the non-central sequence is a forward or reverse chain by comparing these two values ​​(relative to the central sequence, a non-central sequence is a forward chain if it is in the same direction as the central sequence, and a non-central sequence is a reverse chain if it is in the opposite direction; if a non-central sequence has more Kmers in its Kmer dictionary with the forward chain than with its Kmer dictionary with the reverse chain, it is determined to be a forward chain, otherwise it is determined to be a reverse chain).

[0119] When predicting very long sequences, because long sequences contain a large number of Kmers, the probability of prediction errors is extremely low even if only the first half of the sequence is searched in both dictionaries. Therefore, only the first half is searched to accelerate the prediction process. For short sequences, however, because they contain fewer Kmers, less effective information, and are more subject to chance, the probability of prediction errors is higher than for long sequences. Therefore, the entire sequence is traversed to count the number of Kmers identical to those in both directions of the central sequence. Furthermore, since short sequences are typically gene-level sequence data, and their direction is generally the same as the central sequence, this scheme adds a certain weight to the counting of Kmers in the same direction, favoring those identical to the central sequence.

[0120] One-way search for common substrings

[0121] For each non-center sequence, this scheme sequentially uses FM-index to search the BWT data structure of the center sequence to find all common substrings matching the non-center sequence. Initially, due to the short search distance, many common substrings can be found. As the step size gradually increases, it continues until no common substrings are found, then backtracks one step and stops. In this process, it can be determined that the sequence starts from the starting position B_index, and after a step size len, x (x>0) common substrings are found in the center sequence. At the same time, the starting positions of these common substrings in the center sequence (i.e., the values ​​of the corresponding suffix array SA) A_index_1, A_index_2…A_index_x can also be obtained. Next, it is determined whether the length of len is greater than the common substring length threshold threshold_BWT. If so, this information is recorded as common substring information. The recorded common substring information is either a one-to-one or one-to-many pattern depending on the search results. A one-to-many pattern takes the form [B_index, A_index_1, len], [B_index, A_index_2, len]...[B_index, A_index_x]. Then, the search continues to the next step, starting from the position B_index + len. If the length len does not meet the threshold_BWT, the search restarts from the position B_index + 1. This process continues until the non-center sequence reaches its end, at which point the search for the common substring in that direction is complete. The threshold_BWT is pre-set in each iteration and decreases with each iteration. A larger threshold results in fewer x values, and a smaller threshold results in more x values.

[0122] Bidirectional search for common substrings

[0123] For FASTA format output, only one forward common substring search as described above is performed. For MAF (multiple alignment format) format output, this scheme searches in both forward and reverse directions to find possible inverted structure variations. Unlike the forward search, when a reverse search fails, the next search does not start from B_index+1, but from B_index+100. Inverted structure variations are not common in the entire sequence, so the probability of finding a common substring reaching threshold_BWT in the reverse search is low, leading to a higher likelihood of search failure. Single-step movement would be very slow. Increasing the search speed by 100 positions after each failure allows for quick skipping of regions with extremely low probability, improving search efficiency.

[0124] S4. Filtering the longest common substring combination and finding the SV chain.

[0125] The found common substrings will have very high coverage, small positional offsets, and overlaps. Moreover, they contain structural variations such as large-scale copying / loss, rearrangement, and inversions caused by racial differences and diseases like cancer. Therefore, finding the main chain and the structural variation chains is of paramount importance.

[0126] Find the main chain

[0127] For small-scale genomic data, since the number of common substrings found above is limited, the most accurate but time-consuming dynamic programming algorithm can be used to find the main chain when applying this method. For large-scale genomic data, the number of common substrings is enormous, and the dynamic programming algorithm with a time complexity of O(n²) is extremely slow. Therefore, this approach uses the following steps to filter them:

[0128] 1. First, all the retrieved one-to-many homologous intervals are filtered into one-to-one patterns according to the principle of closest relative position, so as to minimize the offset between the two sequences and make the subsequent comparison more uniform.

[0129] 2. Next, we address the overlap problem by treating the common substrings and their common origin intervals as individual intervals. After sorting these intervals by their right ends, we use a greedy algorithm to select the longest non-overlapping sub-interval.

[0130] 3. Finally, the longest increasing subsequence (LIS) algorithm is used to solve the cross-permutation problem, so as to obtain the final combination of the same source intervals.

[0131] like Figure 2 Figure B illustrates an example of finding the main chain within a common substring. The figure shows five homologous intervals (X1-X5) retrieved from the sample sequence (non-center sequence) relative to the reference sequence (center sequence). X1 corresponds to Y1 and Y2, X2 to Y3, X3 to Y8, X4 to Y4, Y5, and Y6, and X5 to Y7. Since X1 and X4 are a one-to-many homologous interval, we first filter X1 and X4 within multiple homologous regions based on their relative positions, making it a one-to-one pattern: X1 selects Y1, and X4 selects Y5. Then, we address the overlap issue, finding that X4(Y5) and X5(Y7) overlap. Since X4 has the longest overall interval length, we discard the homologous region represented by X5(Y7). Finally, the problem of cross-permutation was solved. It was found that X3(Y8) and X4(Y5) cross each other. Since X4 has the longest interval length, the homologous region represented by X3(Y8) was discarded. Thus, the main chain composed of X1(Y1), X2(Y3), and X4(Y5) was obtained.

[0132] In the MAF mode of this scheme, finding common substrings needs to be performed in both the forward and reverse directions. After obtaining the common substrings in both directions, the search for the main chain also needs to be performed in both directions, thus obtaining the longest common substring combination in both directions. Among them, the longest common substring combination in the same direction as the center sequence is the longest and contains the most complete information, while the longest common substring combination in the reverse direction is actually formed by a few inverted structure mutations. The longest common substring combination in the reverse direction is divided into multiple inverted clusters, each corresponding to an inverted structure mutation. These inverted clusters are merged with the forward main chain. In addition, when overlap occurs, this scheme prunes the corresponding inverted clusters, finally obtaining the true global main chain. Once the inverted clusters are integrated into the global main chain, the forward non-homologous regions corresponding to the inverted clusters do not need to be dynamically programmed or divide-and-conquer compared, thereby reducing the overhead of comparing low-quality regions.

[0133] like Figure 2 The diagram in Figure C illustrates the process of merging the forward and reverse main chains. The first step is to locate the forward main chain, indicated by the solid arrow; the second step is to locate the reverse main chain, indicated by the dashed arrow; the third step is to map the reverse main chain to the forward main chain; and the fourth step is to merge the forward and reverse main chains. The solid and dashed arrows in the diagram form the final main chain, and the intervals within it no longer need to be compared; only the gaps require further dynamic planning and comparison.

[0134] Finding SV chains

[0135] Besides the inverted clusters in the reverse main chain, which constitute an SV chain, this scheme uses the common substrings filtered out for finding forward and inverse main chains as the data source for finding other SV chains. This data is used to perform two separate SV chain searches, identifying both forward structural variations and structural variations combined with inverted clusters. Preferably, this scheme uses several thresholds to limit the common substrings in suspected chains, in order to find SV chains composed of these compact, parallel common substrings. Figure 2In D, the specific restriction rules are: Alen>thresh_SV; Blen>thresh_SV; disA<thresh_dis; disB<thresh_dis; all_len / min(ALen,BLen)>thresh_RATIO; abs(disA-disB) / min(disA,disB)<thresh_sub_ratio; abs(ALen-BLen) / min(ALen,BLen)<thresh_SUB_RATIO, wherein thresh_SV is the minimum length of structural variation specified by the user, thresh_dis is the maximum distance between adjacent common substrings, thresh_RATIO is the minimum ratio of the total length of common substrings to the chain length, thresh_SUB_RATIO is the maximum relative error of the total lengths of two sequences in the chain, and thresh_sub_ratio is the maximum relative error of the distance between two adjacent common substrings in the chain.

[0136] After obtaining some compact, parallel and regular structural variation chains, if there is coverage overlap, these chains are further screened. For dual genome alignment, in the alignment result file, as for whether a local alignment result (block) is retained, it is determined whether both the ref region and the sample region in the block are covered by other blocks; if yes, the block is not retained, otherwise the block is retained. Finally, Kband dynamic programming is used to perform alignment filling on the regions between common substrings in the chain to obtain the complete aligned structural variation alignment, which is written into a maf file in the form of blocks.

[0137] S5. BWT recursive divide-and-conquer for non-homologous intervals

[0138] BWT / FM-index can quickly perform conversion compression and common substring retrieval on a central sequence, and the optimal common substring combination matching the central sequence can be screened out through greedy algorithm and LIS problem solving. Therefore, the above method forms a complete chain for solving the optimal common substring combination between two sequences: BWT-FM-LIS.

[0139] For large genome sequences like those of humans, this approach employs a BWT-FM-LIS cyclic divide-and-conquer method. By progressively lowering the common substring length threshold during the retrieval phase (default global threshold_BWT = 15, local cyclic divide-and-conquer threshold_bwt = 5), the flexibility of BWT-FM-LIS is enhanced, allowing it to be applied to finding optimal common substring combinations between sequences of varying lengths and similarities. For extremely long genome sequences (such as human chromosome 1), a longer threshold can be used to quickly split the 248Mb sequence, filtering out the main strand and obtaining shorter but less similar differential substring pairs. These differential substring pairs are then further decomposed using a shorter common substring length threshold. This iterative process continues until the differential substring pairs are segmented into fragments of suitable length for dynamic programming (default threshold_DC = 10000). Finally, a K-band dynamic programming algorithm is used for fine-grained alignment.

[0140] S6. K-band dynamic programming for non-originating intervals

[0141] After the BWT-FM-LIS cyclic divide-and-conquer approach, the alignment of non-homologous regions in the genome is simplified to dynamic programming among many short sequences. This scheme optimizes the dynamic programming from four aspects: scoring rules, constraints, data storage, and underlying system memory characteristics.

[0142] Affine penalty improves comparison quality

[0143] Among various genomic variations, the frequency of occurrence is typically substitution >> deletion >> insertion >> large-segment structural variations. After the initial insertion or deletion, subsequent insertions or deletions are relatively easy to occur. To make the alignment results more closely resemble natural patterns, an affine penalty is used to distinguish between initial and subsequent insertions: the penalty for initial insertions is greater than that for subsequent insertions, ensuring that insertions located close to each other in the alignment results are as adjacent as possible. This helps to obtain more accurate variations and their biological significance, such as... Figure 2 In the diagram, E and -σ represent the initial gap penalty, -e represent subsequent (non-initial gap) penalties, and +s represents the normal alignment score determined by the scoring matrix. Generally, the initial gap penalty is heavier, while subsequent gap penalties are lighter (gap penalties refer to the comparison between A / C / G / T and -, where - represents the gap). The alignment score is determined solely by the scoring matrix: 4[ACGT]*4[ACGT].

[0144] Limiting the scope of dynamic programming saves time and space costs.

[0145] For highly similar sequences, the time and space complexity is reduced from O(n²) to linear O(kn) by limiting the range of dynamic programming, where n represents the length of the sequence. This scheme uses the K-band method to restrict the dynamic programming to the vicinity of the diagonal. It primarily utilizes the fact that the backtracking path for optimal alignment between two highly similar sequences is generally located near the diagonal. Therefore, it is unnecessary to fill and calculate the entire matrix; only the region near the diagonal needs to be calculated, which is called the K-band. The baseline of the band is determined by the difference in length between the two sequences, which is the region between two parallel lines close to the diagonal. Then, using this band as a base, the algorithm expands to the upper right and lower left, starting from 1 and increasing in powers of 2, until the score of each subsequent change is no higher than the score of the previous change, at which point the current dynamic programming ends. Figure 2 F is a specific example of K-band dynamic programming. The length difference between the two sequences is 1, and the width of kband is 1 + 2 * k + 1. Therefore, the light gray band in the middle (located between the two dark gray bands) has k = 0, and the light gray + dark gray band region has k = 1. The maximum score and backtracking path of k do not change from 0 to 1, so the expansion stops when k = 1.

[0146] Path storage simplifies data structure

[0147] The process of dynamic programming involves first calculating the state table, then backtracking to find the optimal path based on the state table. However, adding affine penalty points requires three state tables to store the dynamic programming scores. This solution uses a separate path information table to store the path information during calculation. During backtracking, the optimal path can be found directly from this path information table. Therefore, it's unnecessary to store the scores of the entire state table during calculation; only the current three states need to be stored. Calculations are performed row by row, overwriting the previous row, and the path direction for each state transition is recorded throughout the process. Previously, the three affine penalty state tables were stored using ints, with each int occupying 4 bytes, totaling 12 bytes for three states. The path information table, however, only requires 1 byte for every three states, using only six bits to record the path transition direction for the current state (see...). Figure 2 (G), which means that the memory requirement here will be reduced by about ten times.

[0148] Opening up global space to replace repeated applications and releases

[0149] The Kband algorithm is already at a low level in the program and will be called frequently. Within each Kband method, the dynamic programming matrix requires frequent allocation and deallocation as k increases. Research into the underlying operating system's memory management revealed that frequent memory allocation and deallocation reduces memory operation efficiency, causing memory to accumulate even after being deleted / freed, resulting in a false increase in memory usage. Therefore, this solution allocates the maximum potentially usable memory at once, replacing repeated allocation and deallocation with repeated read / write operations. It abandons the kn space complexity of Kband and chooses n² instead, using the kn time complexity of Kband to optimize performance, thus reducing peak memory usage to an ideal level.

[0150] S7. Data Post-processing

[0151] Star Comparison Integration

[0152] In large-scale multiple sequence alignment of the human genome, using the human reference genome as the central sequence aligns perfectly with the idea of ​​star alignment. Each sequence only needs to be aligned with the reference genome once, and then the multiple sequence alignment is completed by integrating the gap information from each pair of double sequence alignments according to the principle of "once empty, always empty". This is because only the gap information (index 5, inserting 3 gaps) is actually recorded as <5,3>. Since the gaps are few, this recording method does not consume much memory, unlike the entire sequence ACGT---ACGT, which is very memory-intensive. This process updates the original double sequence alignment gap information with some <> that can make the sequences aligned in the multiple sequence alignment. Finally, the storage method is the same as N, hence the following module, which inserts N and gaps back into the actual sequence.

[0153] Integration of degenerate base N

[0154] This scheme has already obtained the insertion information for each sequence through star alignment, but the current alignment is based on the sequence after preprocessing to remove degenerate bases (hereinafter, N represents all degenerate bases). Therefore, the N removed during preprocessing must be integrated into the final alignment. There are three types of indices corresponding to the following three types of strings:

[0155] Original string: The original string with N;

[0156] Input string: A string containing only ACGT characters after removing N;

[0157] Align string: The string after adding gaps to the input string for alignment (alignment is achieved through star comparison).

[0158] 1. The main body of the program does not contain N. When reading in the initial data, N is removed, but it is recorded as a binary vector of (index, num). Later, the information of gap is also obtained as a binary vector of (index, num). The index of N and gap is relative to the "input string".

[0159] 2. Integrate the two insertion vectors into triplet vectors (index, n_num, g_num), and insert them in reverse order according to the "input string" index (because inserting at the end first will not affect the index at the beginning, while inserting at the beginning first will change the index of the subsequent insertions, causing the insertion to be disordered and incorrect). However, ensure that g_num is the upper limit. That is, if there are 4 gaps and 6 N at the same position, only 4 N can be inserted. The remaining two N will be considered later, because the alignment of the "Align string" must remain unchanged.

[0160] 3. For the extra N and the gap in other strings caused by the extra N in a certain string (the principle of same-star comparison), insert them in reverse order according to the index of the "Align string".

[0161] 4. Finally, special cases were handled. If three sequences were aligned, and a certain position initially had no gap, but two strings had N, then gaps would inevitably need to be inserted at that position in all the other strings except the original sequence. This would result in an extra column consisting entirely of gaps. Later, the insertion vector for the "Align string" was modified to avoid this unnecessary insertion. For example... Figure 2 In the example, H represents a specific instance of deleting N and inserting N in a sequence of 3. Inserting N back into the first sequence results in two light gray gaps ("-") in the remaining sequences. Inserting N back into the second sequence results in two dark gaps in the remaining sequences. In the end, a column is composed entirely of gaps, which is unreasonable and needs to be deleted.

[0162] Insert Module

[0163] In this scheme, the insertion operation on the original sequence is not performed directly. Instead, the insertion operation is transformed into an operation on the pair vectors of insertion N and insertion gap. This design postpones the time-consuming actual insertion operation to the end, avoiding wasted time due to unsuccessful alignment after the actual insertion for some reason.

[0164] While linked lists are the most efficient data structure for insertion operations, it's generally difficult to allocate such a large amount of memory to the operating system for such a large sequence and the insertion of information. Furthermore, linked lists have low space utilization. In 32-bit compilation mode, each node of a linked list consists of a 1-byte unsigned char data block and a 4-byte pointer; while in 64-bit compilation mode, the pointer is 8 bytes. Therefore, the memory utilization of linked lists is only 20% or 11%. In addition, linked lists require frequent memory allocation and deallocation operations, leading to reduced efficiency of memory operations by the underlying operating system.

[0165] To address this issue, this solution employs a different insertion approach: first, the length after insertion is calculated; then, a vector of appropriate size is created, and the original string information and the insertion information are written into this vector sequentially, thus implementing the insertion operation. In this way, the time complexity of the insertion operation for each sequence is only O(n) of a linear traversal.

[0166] S8. Structural Variation and the Merging of MSA

[0167] FSATA is the input format and one of the output formats in this solution. MAF format is another output format in this solution. Compared with FASTA format, it can store alignments of structural variations such as inversions, transpositions, and repetitions. This solution optimizes the functionality of multiz for merging two MAF files, enabling it to be used for merging MAF and FASTA files.

[0168] Multiz itself does not have a sorting function, and it reads the entire file into memory for linked list merging each time. In contrast, this solution sorts the MAF file before merging, saves the sorting results, reads each block one by one and performs the multiz merging operation, and writes out the merged result. This avoids keeping the entire MAF file resident in memory and is more suitable for processing large-scale genomic data.

[0169] Furthermore, multiz can only process the merging of two files in a single thread. In contrast, this solution adopts a layer-by-layer merging method, where the merging within each layer is independent of each other. Therefore, multi-threaded parallel processing can be used, which is suitable for the solution to process multiple MAF files for multiple sequence / genome alignment.

[0170] For the final step of merging structural variation MAF files and multiple sequence alignment result FASTA files, this approach developed an integration module with no information loss that is far more efficient than multiz (e.g., Figure 2(I) First, the blocks in the MAF of SV are sorted according to the center sequence index, and sequences in the blocks with alignment quality exceeding those in the original FASTA are selected. Then, the FASTA alignment matrix is ​​column-cut based on the start and end points of each sequence in the block. Finally, the blocks of structural variants in SV and the corresponding blocks after FASTA cutting are integrated based on star alignments. At the same time, the alignments of structural variants in the original FASTA are deleted. This gives the final MAF format multi-genome alignment result that integrates structural variants.

[0171] To demonstrate the advantages of this solution over existing technologies, a comparative analysis was conducted between this solution and some existing technologies. The software based on this solution is referred to as Halign-G.

[0172] Figure 3 The diagram shows a comparative analysis of Cactus, Parsnp, and Halign-G on 24 chromosome datasets from 18 individuals. Figure A represents the distribution of the M-scores for the three software programs. The upper subplot shows an expanded view of the 0.9-1 range. The three box lines within each vertical frame represent the distributions of Cactus, Parsnp, and Halign-G, respectively. For clarity, 'a' represents the distribution of Halign-G, 'b' represents Parsnp, and 'c' represents Cactus. It can be seen that Parsnp and Cactus have no results on datasets 13, 14, 15, and 22. Cactus's M-score surpasses Halign-G on datasets 5, 7, 8, 16, and 17, while Halign-G has the highest M-score on the remaining 15 datasets. B shows the time and memory consumption of the three software programs in the experiment. The height of the bubble represents the time duration, and the size of the bubble represents the memory usage. Halign-G is the darkest, followed by Cactus, and Parsnp is the lightest. Halign-G is always on the far left of each row, indicating the shortest time consumption and the smallest circle with the least memory usage. It can be seen that Halign-G has the lowest time and space cost.

[0173] Figure 4The comparison analysis shows seven software programs—Cactus, Parsnp, Mugsy, Kalign, MLAGAN, PMauve, and Halign-3—against Halign-G's FSATA and MAF modes across five datasets. A represents the distribution of the M-scores for the nine software programs. B shows the time and memory consumption of the nine software programs in the experiments; the height of the bubble represents the time duration, and the size of the bubble represents the memory usage. For the Mitochondrial-genome, Streptococcus-pneumoniae, and Escherichia-coil datasets, Halign-G's MAF mode performed best in terms of average M-score. Notably, the FASTA mode, with relatively lower alignment quality, also outperformed Halign-3 in terms of M-score at the multiple sequence alignment level. Furthermore, Halign-G had the lowest memory usage among all software programs across all five datasets. The distribution and bubbles in the figure were originally distinguished by color to differentiate the results of each software. However, to avoid the grayscale from affecting the understanding, x is now used to represent the MAF result of Halign-G, y is used to represent the FSATA result of Halign-G, and the rest are the results of the control group, which have little impact on the understanding, so they are not marked separately.

[0174] Figure 5 The image shows a comparison of structural variations identified by Cactus, Parsnp, Mugsy, Pmauve, and Halign-G across nine simulated datasets. Each pentagon in the image represents a standard answer site for a structural variation, divided into five small triangles based on the five software programs. A lit triangle indicates that the corresponding software program has identified that site. In the pentagons, Cactus corresponds to the upper right triangle, Parsnp to the lower right triangle, Halign-G to the bottom triangle, Mugsy to the upper left triangle, and Pmauve to the lower left triangle.

[0175] Experimental results show that among a total of 108 structural variation sites, Pmauve found no sites, Cactus found 21 sites, Parsnp found 35 sites, Mugsy found 52 sites, while Halign-G found 107 sites, missing only one site in the seventh set of simulation data. This fully demonstrates Halign-G's powerful ability to identify structural variations.

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

Claims

1. A multi-genome alignment method applicable to large-scale genomes, characterized in that, The method includes: S1. Perform data preprocessing on the genome sequences to be aligned and determine the center sequence; S2. Establish a BWT index for the center sequence to obtain the BWT data structure, and build a prediction dictionary; S3. Use a prediction dictionary to predict the direction of each non-center sequence; S4. Based on the common substring length threshold, perform a forward common substring retrieval for each non-center sequence in the BWT data structure of the center sequence; S5. Filter the longest common substring combination as the main chain of the corresponding non-center sequence. The common substring information is either a one-to-one mode or a one-to-many mode according to the search results. The methods for finding the main chain include: S51. First, all common substrings retrieved for the corresponding non-center sequence are filtered into a one-to-one mode according to the principle of closest relative position; S52. Find overlapping common substrings and intersecting common substrings. For each pair of overlapping common substrings and each pair of intersecting common substrings, retain the longest corresponding homologous interval in the non-center sequence and discard the other corresponding homologous intervals in the non-center sequence to filter out the longest homologous interval combination. S6. Filter the main chain and determine whether the length of the different substrings meets the requirements of dynamic programming. If yes, continue to S7; otherwise, reduce the threshold for the length of the common substring and repeat S4-S6 to further decompose the different substrings. S7. Use the Kband dynamic programming method to perform a fine comparison of the difference substring pairs: First, initialize the state table; Then, the score of the current row is calculated using affine penalty, and the path transition direction is filled into the state table based on the relationship between the scores of the current row and the previous row. Finally, the optimal path is traced back from the state table; In Kband dynamic programming, Affine penalty is used to differentiate between first and second insertions, so that the penalty for the first insertion is greater than that for the second insertion; Save time and space costs by limiting the scope of dynamic programming: restrict dynamic programming to the vicinity of the diagonal, and only calculate the region near the diagonal; the baseline of the strip is determined by the difference in length between the two sequences, and based on the strip, expand to the upper right and lower left, starting from 1 and increasing in powers of 2, until the score of the next change is no higher than the score of the previous change. Simplify the data structure by storing path information: During the calculation process, a separate information table is set up to store the path information. When backtracking, the best path is found directly by looking up the path information table. During the calculation process, the current three states are stored, and then each row is calculated to overwrite the previous row. The path transfer direction of each state is recorded throughout the process. Instead of repeatedly allocating and freeing memory, allocate global space: allocate the maximum memory that may be used in one go, and replace repeated memory allocation and freeing with repeated read and write operations; S8. Output the results of multiple sequence alignment.

2. The multi-genome alignment method applicable to large-scale genomes according to claim 1, characterized in that, In step S1, Data preprocessing includes: recording the positions of degenerate bases in binary form and removing them along with the gaps; The center sequence can be determined by either selecting the longest sequence after removing degenerate bases, or by specifying it by the user. In step S2, two hash dictionaries are built for the forward and reverse chains of the center sequence to obtain the prediction dictionary. The prediction dictionary includes the forward chain Kmer dictionary and the reverse chain Kmer dictionary of the center sequence. Step S2 also includes storing the central sequence using a binary compression storage method; In step S3, the number of Kmers shared by each non-central sequence is found in the forward chain Kmer dictionary and the reverse chain Kmer dictionary of the central sequence, and the direction of the non-central sequence is determined accordingly.

3. The multi-genome alignment method applicable to large-scale genomes according to claim 1, characterized in that, In step S4, the forward common substring retrieval method is as follows: S1. Starting from the beginning position B_index of the non-center sequence, the initial len is 0. Continuously increment len ​​by 1 to search for all common substrings in the center sequence until there are no common substrings. Backtrack one step to get x common substrings under the current step length len. S2. Determine whether len is greater than or equal to the common substring length threshold. If so, record the common substring information found under that len, and set B_index=B_index+len, then return to step S1 to repeat the search. If not, then set B_index = B_index + 1 and return to step S1 to repeat the search; S3. The above process continues to retrieve the sequence until the end of the non-center sequence.

4. The multi-genome alignment method applicable to large-scale genomes according to claim 2, characterized in that, In step S7, the difference substring pairs are finely compared by star comparison to obtain the insertion information of each sequence; Furthermore, in step S7 of this method, post-processing is performed on the results after fine comparison: The degenerate bases removed during preprocessing are integrated into the final alignment; The degenerate base insertion process is as follows: first, the length after insertion is calculated; then, a vector of the corresponding size is allocated; and the original string information and the insertion information are written into the vector in sequence to realize the insertion operation.

5. The multi-genome alignment method applicable to large-scale genomes according to any one of claims 1-4, characterized in that, When step S8 is required to output the structural variation alignment results, step S4 performs a bidirectional common substring search on each non-center sequence to find possible inverted structural variations. Step S5 will obtain a forward main chain composed of the longest common substrings in the forward direction and a reverse main chain composed of the longest common substrings in the reverse direction. The forward main chain and the reverse main chain will be merged to obtain the final main chain of the corresponding non-center sequence. The merging process includes: Find the forward main chain and the reverse main chain; map the reverse main chain to the forward main chain; merge the forward main chain and the mapped reverse main chain.

6. The multi-genome alignment method for large-scale genomes according to claim 5, characterized in that, In bidirectional retrieval, the forward retrieval method is: S1. Starting from the beginning position B_index of the non-center sequence, the initial len is 0. Continuously increment len ​​by 1 to search for all common substrings in the center sequence until there are no common substrings. Backtrack one step to get x common substrings under the current step length len. S2. Determine whether len is greater than or equal to the common substring length threshold. If so, record the common substring information found under that len, and set B_index=B_index+len, then return to step S1 to repeat the search. If not, then set B_index = B_index + 1 and return to step S1 to repeat the search; S3. The above process continues until the end of the non-center sequence; The reverse retrieval method is: S1. Starting from the beginning position B_index of the non-center sequence, the initial len is 0. Continuously increment len ​​by 1 to search for all common substrings in the center sequence until there are no common substrings. Backtrack one step to get x common substrings under the current step length len. S2. Determine whether len is greater than or equal to the common substring length threshold. If so, record the common substring information found under that len, and set B_index=B_index+len, then return to step S1 to repeat the search. If not, then set B_index = B_index + 100 and return to step S1 to repeat the search; S3. The above process continues to retrieve the sequence until the end of the non-center sequence.

7. The multi-genome alignment method for large-scale genomes according to claim 6, characterized in that, This method also includes: The common substrings filtered out by the forward main chain and the reverse main chain are used as the data source for finding other structural variant chains; Two search operations were performed on the data source to find the positive structural variation chain and the structural variation chain combined with the inversion; Selective retention processing is applied to overlapping structural variant chains; In step S7, the dynamic programming method is used to compare and fill the regions between common substrings of the structural mutation chain to obtain the entire aligned structural mutation alignment, and write it to the maf file; The results of multiple sequence alignments are written to a FASTA file; In step S8, the maf file and fasta file are merged to output the final multi-genome alignment results: First, sort the blocks in the MAF file according to the center sequence index and then filter the sequences in the block whose alignment quality exceeds that of the original FASTA file. Then, based on the start and end points of each sequence in the block as cutting points, the FASTA alignment matrix is ​​column-cut. Finally, based on the star comparison, the blocks of the MAF file and the corresponding blocks after the FASTA file are cut are integrated, and the structural variations in the original FASTA file are deleted. Before merging, the MAF files are sorted and the sorting results are saved. Then, each block is read in and merged, and the merged result is written out. The merging process adopts a layer-by-layer merging method to ensure that the merging within each level is independent of each other.

Images