A fast sequence alignment method driven by rare matching anchors for highly repetitive genomic regions

Abstract

The application provides a rare matching anchor point driven fast sequence alignment method for highly repetitive genomic regions, which solves the problems of insufficient alignment accuracy, excessive running resource consumption and difficulty in forming biologically credible insertion and deletion statistics for highly repetitive regions such as centromeres, and uses rare matching fragments as anchors, combines index structure, segmented recursion and efficient and accurate alignment, realizes low memory and scalable double sequence global alignment, and outputs alignment results that can be used for insertion and deletion statistics, repeat expansion analysis and structural difference analysis. The application has the advantages of high alignment quality, good scalability and the like.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and comparative genomics, specifically relating to a method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions. Background Technology

[0002] With the rapid accumulation of high-quality long-read assembly and pan-genome data, researchers increasingly need to make precise comparisons of centromeres, telomeres and other highly repetitive regions in different individuals in order to characterize the amplification, deletion and complex rearrangement patterns of repetitive arrays.

[0003] However, sequences in highly repetitive regions often contain a large number of approximately repetitive units. The common alignment strategy of "short fragment seed localization + extension" is prone to generating a large number of ambiguous matches, leading to unstable alignment chains, blurred structural difference boundaries, and even difficulty in providing continuous and complete global alignment. On the other hand, traditional precise global alignment based on full matrix dynamic programming has extremely high time and space overhead for millions of base pairs, making it difficult to meet the requirements of engineering feasibility. Although existing methods can achieve good results in non-repetitive regions, they still suffer from insufficient alignment accuracy, excessive runtime resource consumption, and difficulty in generating biologically reliable insertion and deletion statistics in highly repetitive regions such as centromeres.

[0004] Therefore, there is an urgent need for a fast and accurate alignment scheme for highly repetitive and ultra-long sequences: one that can achieve global alignment with controllable memory usage and preserve as much of the true evolutionary event signal as possible. Summary of the Invention

[0005] The purpose of this invention is to address the above-mentioned problems by providing a rare matching anchor-driven fast sequence alignment method for highly repetitive genomic regions. This method utilizes rare matching fragments as anchors, combines index structures, segmented recursion, and efficient and accurate alignment to achieve low-memory, scalable global alignment of two sequences, and outputs alignment results that can be used for insertion / deletion statistics, duplication / amplification analysis, and structural difference resolution.

[0006] To achieve the above objectives, the present invention employs the following technical solution: a rare matching anchor-driven fast sequence alignment method for highly repetitive genomic regions, comprising the following steps: Step 1: Input and Preprocessing. Obtain the first sequence S1 and the second sequence S2 to be aligned, and standardize S1 and S2. The input consists of two nucleic acid sequences, which can come from a single sequence file or be obtained by splicing multiple sequence fragments. Before entering the anchor point search, the system performs uniform preprocessing on the input, removing non-base characters that may interfere with index construction to ensure that the output alignment can trace back to the original sequence coordinates. The input FASTA / plain text sequence is scanned character by character: (1) Convert to uppercase; (2) For characters other than {A,C,G,T,N}: discard if they are separators / newlines / spaces; map them to N if they are IUPAC fuzzy bases. The final output is the globally aligned result of the double sequence in FASTA format.

[0007] Step 2: Rare anchor point discovery. Construct an index structure containing the spliced ​​sequence S1 and S2. Based on the index structure, identify one or more common segments that appear at least once in both S1 and S2 and whose total number of occurrences in the spliced ​​sequence meets a preset rare threshold. The occurrence position of each common segment in S1 and S2 constitutes a candidate anchor point. Since rare matching segments are used as anchor point candidates, to efficiently discover such segments, the system constructs an index and auxiliary array suitable for common segment retrieval for the two sequences, and calculates the suffix sorting and the longest common prefix information on the concatenated string. By performing interval scanning on the longest common prefix information, the system can find common segments that appear simultaneously in both sequences, occur infrequently, and reach a certain length, and uses these as rare anchor point candidates. Unlike strategies that rely on common short seeds, rare anchor points emphasize discriminative power and can provide more stable localization criteria in highly repetitive backgrounds.

[0008] Step 3: Anchor chain screening and recursive segmentation. For all candidate anchors obtained in Step 2, chain screening is performed based on collinearity constraints and a predefined scoring function to obtain a collinear anchor main chain with the highest score. Using the anchors in the anchor main chain as boundaries, S1 and S2 are segmented into multiple independent sub-interval pairs. The alignment framework using rare matching fragments as core anchors is based on the idea that in highly repetitive regions, widely occurring short fragments often fail to provide effective localization, while near sequence differences, mutations, or array boundaries, there are still less frequent and more discriminative common fragments. Using these rare common fragments as anchors, more reliable alignment relationships can be established in repetitive backgrounds. Based on this idea, rare anchors are first identified between two sequences to be aligned, and a collinear, consistent anchor chain is obtained through chain-like screening. Then, using the anchor chain as a framework, the alignment task is globally divided into multiple relatively independent sub-regions.

[0009] Step 4: Recursive processing. For each sub-interval pair, repeat steps S2 and S3 until the recursion termination condition is met, resulting in multiple terminating sub-intervals. The anchor point discovery and filtering process is recursively repeated for each sub-interval until the sub-interval length falls within a range that can be precisely processed. This divide-and-conquer strategy decomposes the global alignment of ultra-long sequences into multiple sub-problems of manageable size, enabling the algorithm to maintain accuracy while possessing better scalability and lower memory pressure.

[0010] Step 5: Precise global comparison and completion. For each terminating sub-interval and the gap segment between anchor points, perform precise global comparison based on the wavefront advancement algorithm. The wavefront advancement algorithm uses a segmented gap cost model for scoring. After the anchor point main chain is determined, the system performs precise global alignment on the gaps between anchor points and the edge segments not yet covered by anchor points to generate continuous and complete alignment results. To efficiently complete accurate sequence alignment within the sub-intervals divided by anchor points, a wavefront-progressive alignment algorithm is adopted.

[0011] In this algorithm, the alignment process does not involve a traditional exhaustive traversal of the entire alignment matrix. Instead, it proceeds layer by layer according to the priority of the current alignment score, and dynamically maintains a set of alignment states that are most likely to form the optimal path. Each state represents the current position on the reference sequence and the query sequence, as well as the corresponding alignment decision direction.

[0012] During the process, the algorithm starts from the initial position and generates candidate alignment paths according to three types of operations: matching / mismatching, insertion, and deletion, while updating possible advancement fronts in real time. For each score level or stage, the algorithm prioritizes expanding states closer to the main diagonal direction from the wavefront, while pruning states that deviate too far, thus avoiding unnecessary computational overhead. This strategy of advancing by score and expanding by wavefront allows the algorithm to concentrate resources on exploring regions that are more likely to belong to the globally optimal alignment path, without wasting time on obviously unreasonable matching decisions.

[0013] When a wave of advancement reaches the termination position of the reference sequence and the query sequence, the algorithm stops advancing and backtracks along the recorded state source information to generate the final alignment path and operation sequence. In this way, the alignment strategy not only guarantees the globally optimal and accurate alignment result, but also significantly reduces the time and space overhead required by traditional dynamic programming.

[0014] To further enhance the biological rationality of the expression of insertion and deletion events, a segmented gap cost concept is introduced in the precise alignment stage. That is, different growth laws are used for the cost models of missing gaps and long gaps, so that the algorithm can give an alignment path that is more in line with the actual evolution process when facing events such as large segment deletions and large segment duplications that are common in repeated arrays.

[0015] Step 6: Output the results. Combine the alignment results of all segments in Step 5 with the anchor point information in Step 3 to generate and output the complete global alignment results of S1 and S2.

[0016] For the task tree formed by recursive partitioning, merge the anchor points and the comparison results of each segment according to preorder traversal to obtain the final alignment result. The output is a continuous and complete double-sequence global alignment result, which can further derive structural difference events and coordinate mapping information, facilitating subsequent comparison of centromeric repeat arrays, estimation of indel length distributions, and analysis of repeat structure evolution for pan-genome research.

[0017] In the above rare matching anchor point-driven fast sequence alignment method for highly repetitive genomic regions, the construction of the index structure and the identification of rare matching fragments in step 2 specifically include: Construct the concatenated sequence S = S1 || $ || S2, where $ is a delimiter that does not belong to {A, C, G, T, N}; Construct the suffix array SA, the longest common prefix array LCP, and the inverse suffix array ISA for the concatenated sequence S; By scanning the LCP-interval of the LCP array, identify the common fragments that appear at least once in both S1 and S2 and whose total number of occurrences in S meets the rare threshold as rare matching fragments.

[0018] Let the two input nucleic acid sequences be the reference sequence S1 and the query sequence S2 respectively. Construct the concatenated sequence S = S1 || $|| S2, where $ is a delimiter and $ does not belong to {A, C, G, T, N}. Define the LCP-interval as the interval [i..j] such that min(LCP[i..j]) = l_min, and the boundaries satisfy LCP[i - 1] < l_min and LCP[j + 1] < l_min. For a certain LCP-interval, [i..j] covers n = j - i + 1 LCP elements, then the corresponding common prefix fragment appears n + 1 times in S.

[0019] In the above rare matching anchor point-driven fast sequence alignment method for highly repetitive genomic regions, the rare matching fragment t in step 2 is defined as a common fragment that simultaneously satisfies the following conditions: Condition 1: c_S1(t) ≥ 1 and c_S2(t) ≥ 1, where c_X(t) represents the number of times the fragment t appears in the sequence X; Condition 2: 2 ≤ c_S(t) ≤ C_max, where c_S(t) is the total number of occurrences of the fragment t in the concatenated sequence S, and C_max is a preset maximum occurrence threshold.

[0020] In the above rare matching anchor - driven fast sequence alignment method for highly repetitive genomic regions, in step 2, the search for rare matching fragments adopts a stepped strategy from strict to loose: first search for common fragments with longer length and fewer occurrences as candidate anchors. When the anchors found so far cannot cover the region to be aligned, gradually lower the length threshold and / or raise the maximum occurrence threshold C_max, and continue the search, taking into account both speed and coverage.

[0021] In the above rare matching anchor - driven fast sequence alignment method for highly repetitive genomic regions, the anchor chain screening and recursive segmentation in steps 3 and 4 specifically include: Step 31: For each rare matching fragment t obtained in step 2, combine any occurrence position x of it in S1 with any occurrence position y in S2 to form an anchor A=(x,y,w), where w = |t|, and the anchor represents the matching correspondence between the S1 interval [x, x + w - 1] and the S2 interval [y, y + w - 1]; Step 32: Define the collinear connectable constraint: Given two anchors Aj=(xj,yj,wj) and Ai=(xi,yi,wi), if xj + wj ≤ xi and yj + wj ≤ yi are satisfied, then Aj→Ai is said to be collinearly connectable; Step 33: Define an anchor scoring function, including a benefit term and a penalty term: The benefit term α(i)=li / min(c_S1(ti), c_S2(ti)), where li is the length of the fragment ti corresponding to the anchor Ai, and c_X(ti) represents the number of occurrences of the fragment ti in the sequence X; The penalty term β(j,i)=0 when |(yi - yj)-(xi - xj)| = 0; otherwise β(j,i)=2·log2(|(yi - yj)-(xi - xj)|); Step 34: Sort all candidate anchors in ascending order of the x - coordinate, and use the dynamic programming recurrence formula f(i)=max_{j < i, Aj→Ai connectable} ( f(j) + max(α(i)-β(j,i), ε) ) for chain screening to obtain the main chain of collinear anchors with the highest score, where ε is a preset positive threshold; Step 35: Using the anchors in the anchor main chain as boundaries, divide S1 and S2 into multiple independent sub - interval pairs, and recursively execute steps 2 to 35 for each sub - interval pair; Step 36: The recursive termination condition is: within the current sub - interval pair S1[l1..r1] and S2[l2..r2], no common fragment that satisfies the definition of the rare matching fragment in step 2 can be found.

[0022] In the above-mentioned rare matching anchor point-driven fast sequence alignment method for highly repetitive genomic regions, in order to more reasonably characterize the evolutionary characteristics of the coexistence of short gaps and long gaps commonly found in regions such as centromeres, a piecewise growth method is adopted for the gap cost: a more sensitive cost growth is adopted for short gaps to suppress excessive fragmentation, and a relatively smooth cost growth is adopted for long gaps to accommodate large fragment insertions and deletions, so that the output alignment has more interpretability in terms of insertion and deletion statistics. The piecewise gap cost model in step 5 is defined as the minimum of two affine penalty functions: g_2p(k) = min( O1 + E1·(k-1), O2 + E2·(k-1) ), where k≥1 is the length of consecutive gaps, O1 and E1 are the gap opening penalty and extension penalty of the first affine penalty, O2 and E2 are the gap opening penalty and extension penalty of the second affine penalty, and O2>O1, E2<E1.

[0023] In the above-mentioned rare matching anchor point-driven fast sequence alignment method for highly repetitive genomic regions, when processing the recursively generated sub-intervals, the global index structure constructed in step 2 is reused to derive the local index of the sub-intervals. The derivation methods include: using the inverse suffix array ISA to extract and sort the global rankings of positions within the sub-interval to obtain the sub-interval suffix array, and calculating the sub-interval LCP array by performing range minimum query RMQ on the original LCP array.

[0024] To avoid reconstructing the index for sub-intervals, a strategy of constructing once and reusing throughout is adopted. For each position p within the sub-interval, take its ranking ISA[p] in the original SA; after sorting these rankings and mapping them back to the original SA, the suffix order of the sub-interval (sub-interval SA') is obtained. Utilize the property LCP(u,v)=min(LCP[u+1..v]). For adjacent suffixes in the sub-interval SA' with subscripts u<v in the original SA, their adjacent LCP' can be obtained by performing RMQ (range minimum query) on the original LCP. RMQ can be implemented with O(1) query through linear preprocessing such as a block sparse table, thus supporting the rapid construction of a large number of sub-interval-derived indexes.

[0025] In the above-mentioned rare matching anchor point-driven fast sequence alignment method for highly repetitive genomic regions, steps 2 to 5 are parallelized to cover two types of work: rare anchor point search and sub-interval exact alignment. A fixed-size thread pool T is adopted, and a work queue is used for scheduling.

[0026] First, two types of tasks are defined: the anchor task AnchorTask, which is used to execute steps 2 and 3; the alignment task AlignTask, which is used to execute step 5; The scheduling strategy prioritizes AnchorTasks to generate more subproblems quickly. When a terminating sub-interval is generated, the corresponding AlignTask is submitted to the thread pool.

[0027] AnchorTask(interval): Derive / reuse indexes on the interval, search for rare matching fragments, generate candidate anchors and perform chaining DP to obtain the main chain; then split the main chain to generate new subtasks.

[0028] AlignTask(interval): Performs WFA(two-piecegap) on the interval that triggers the recursion termination condition or has no anchor point to obtain the CIGAR of that interval.

[0029] Final result merging: The task tree formed by recursive segmentation is traversed in preorder and the anchor points are compared with the results of each segment to obtain the final alignment result.

[0030] In the above-mentioned rare matching anchor-driven fast sequence alignment method for highly repetitive genomic regions, the standardization process in step 1 includes: converting sequence characters to uppercase, discarding non-base separators or whitespace characters, and mapping IUPAC fuzzy bases to N.

[0031] In the above-mentioned rare matching anchor-driven fast sequence alignment method for highly repetitive genomic regions, the complete global alignment result output in step 6 includes the CIGAR string, the list of structural difference events, and the coordinate mapping relationship between S1 and S2, which facilitates subsequent centromere repeat array comparison, insertion and deletion length distribution estimation, and repetitive structure evolution analysis for pan-genome research.

[0032] Compared with the prior art, the present invention has the following advantages: 1. Design the overall process for scenarios with highly repetitive and ultra-long sequences such as centromeres. Use rare matching fragments as highly discriminative anchor points to establish more reliable positioning relationships in repetitive backgrounds, thereby improving alignment quality and reducing mismatches and discontinuities caused by ambiguous matching.

[0033] 2. By using the "anchor main chain + recursive segmentation" organization method, the global alignment of ultra-long sequences is decomposed into multiple sub-problems of controllable size. In the recursive process, indexes and auxiliary information are reused, which significantly reduces the overhead of repeated construction. This allows the algorithm to maintain accuracy while having better scalability and lower memory pressure.

[0034] 3. In the precise completion stage, a wavefront-driven alignment strategy is introduced, combined with a segmented gap cost model, which enables the algorithm to handle large-segment insertion / deletion and replication / amplification events commonly found in highly repetitive regions more naturally, thereby improving the biological rationality of insertion / deletion statistics and structural difference analysis.

[0035] 4. It supports parallel processing of key stages such as anchor point discovery and precise gap comparison, which can make full use of multi-core computing resources to improve overall throughput. At the same time, it controls additional memory overhead through unified scheduling and intermediate structure reuse, making it suitable for batch sample or long fragment alignment tasks.

[0036] 5. The output is a continuous and complete global alignment result of double sequences, and can be further derived with structural difference events and coordinate mapping information, which facilitates subsequent comparison of centromere repeat arrays, estimation of insertion and deletion length distribution, and evolutionary analysis of repeat structures for pan-genome research. Attached Figure Description

[0037] Figure 1 This is a schematic diagram of the overall process of the present invention. Detailed Implementation

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

[0039] like Figure 1 As shown, this embodiment demonstrates the core workflow of a rare matching anchor-driven fast sequence alignment method for highly repetitive genomic regions. Two short nucleic acid sequences are input to verify the effectiveness of the rare anchor identification and recursive segmentation mechanism. The specific workflow includes the following: Input and Preprocessing: The input sequences are S1 = “AGAGATTGAGCGAG” and S2 = “AGCGAGATAGCGAG”. First, preprocessing is performed to convert both sequences to uppercase, check and remove illegal characters, and establish the original coordinate mapping.

[0040] Rare anchor point discovery: Construct the spliced ​​sequence S = “AGAGATTGAGCGAG$AGCGAGATAGCGAG”, and construct suffix array SA, LCP array, and inverse suffix array ISA for S. Common segments are identified by scanning the LCP-interval. “AG” is a high-frequency repeating unit, appearing frequently in S, and does not meet the rare threshold (assuming C_max is set to 3, c_S(“AG”) = 5 > C_max). “AT” is relatively rarer, appearing once in S1 and once in S2, with c_S(“AT”) = 2, satisfying the condition 2 ≤ c_S(t) ≤ C_max. Therefore, “AT” is identified as a rare matching segment, and its occurrence position constitutes a candidate anchor point.

[0041] Anchor chain screening and recursive segmentation: Since there is only one rare anchor in this example, it is directly used as the main anchor chain. Using this anchor as the boundary, S1 and S2 are segmented into a left sub-interval (S1[1..3]“AGA” and S2[1..3]“AGC”) and a right sub-interval (S1[5..14]“GAGCGAG” and S2[5..14]“TAGCGAG”). Rare anchor discovery is recursively performed on the left and right sub-intervals. After detection, there are no common segments in the sub-intervals that meet the rare threshold condition, which meets the recursion termination condition, and the process proceeds to precise comparison.

[0042] Precise global alignment and completion, and result output: The wavefront advancement algorithm is executed on the left and right terminating sub-intervals respectively, and a segmented gap cost model is used for scoring (parameter settings: O1=2, E1=1, O2=5, E2=0.5). The alignment result of the left sub-interval is "AGA" and "AGC" corresponding to CIGAR "2M1X", and the alignment result of the right sub-interval is "GAGCGAG" and "TAGCGAG" corresponding to CIGAR "1X6M". Finally, the anchor point information and the alignment results of each segment are merged to obtain the global CIGAR "2M1X·2M (anchor point 'AT')·1X6M", and the global CIGAR string, the list of structural difference events (mismatch G↔C exists at position 3 of S1, mismatch G↔T exists at position 5 of S2), and the coordinate mapping relationship are output. Example 2

[0043] This embodiment demonstrates the application effect of the present invention in highly repetitive, ultra-long genome regions such as centromeres, and specifically includes the following process: Input and Preprocessing: The input sequences are the centromere region sequence S1 of the human X chromosome and another haplotype sequence S2 from the same region. After preprocessing and standardization, the original coordinate mapping is established.

[0044] Index Construction and Rare Anchor Discovery: A spliced ​​sequence S = S1 + "$" + S2 was constructed. A linear-time algorithm was used to construct the suffix array SA, LCP array, and inverse suffix array ISA, with a memory usage of approximately 5-6 times the length of the original sequence. A stepwise search strategy was adopted: In the first round, the minimum length L_min was set to 50 and C_max to 5. Scanning the LCP-interval identified approximately 1,200 candidate anchors, with lengths ranging from 50 to 200 bp. These were segments that appeared 2-5 times in both sequences, mainly located in the boundary regions of α satellite monoliths or mutation hotspots. For intervals not covered by anchors, the length threshold was gradually reduced to L_min = 30 and 20, and C_max was increased to 10 and 15 for supplementary searches. Finally, approximately 3,500 candidate anchors were obtained.

[0045] Anchor chain selection and recursive segmentation: All candidate anchors are sorted by their x-coordinates. Dynamic programming is used for chain-like selection. The benefit term α(i) = li / min(c_S1(ti), c_S2(ti)) ensures that longer and rarer anchors receive higher benefits. The penalty term β(j,i) = 2·log2(|(yi-yj)-(xi-xj)|) penalizes collinearity deviations. The optimal anchor chain containing approximately 800 anchors is obtained recursively. Using this as the main boundary, the 1.2Mb sequence is segmented into 801 independent sub-interval pairs. Rare anchor discovery and selection are recursively performed on each sub-interval until no rare anchors can be found.

[0046] Sub-interval index reuse: When processing recursive sub-intervals, local indexes are quickly derived using global SA, ISA, and LCP: the ISA values ​​of the positions within the sub-interval are extracted, sorted, and mapped back to SA to obtain the sub-interval SA'. The positions of adjacent suffixes in SA' in the original SA are obtained by querying the original LCP through RMQ to obtain the sub-interval LCP'. This method avoids the overhead of repeatedly constructing suffix arrays.

[0047] Accurate global alignment completion and parallel acceleration: Wavefront-progressive alignment is performed in parallel on approximately 2,500 terminating sub-intervals using a multi-threaded thread pool and work queue scheduling. AnchorTasks are prioritized to quickly generate more sub-problems, and AlignTasks are delivered to the terminating sub-intervals using a segmented gap cost model. The alignment results can accurately identify array differences in α-satellite repeat units, such as expansion from 171bp to 171bp×n, insertion / deletion events at repeat array boundaries, and monosomal sequence differentiation caused by mutations.

[0048] Results and performance evaluation: Merge all sub-interval alignment results and anchor point information to generate a complete global alignment, and output the CIGAR string, a list of structural difference events, and coordinate mapping relationships. Example 3

[0049] This embodiment demonstrates the specific scheduling process of the parallel acceleration framework of the present invention, as follows: Parallel framework initialization: Create a fixed-size thread pool T=8 and a work queue, support task priority scheduling, and define two task types: AnchorTask (high priority) is used to perform rare anchor discovery and chained filtering, and AlignTask (low priority) is used to perform precise global comparison.

[0050] Initial task delivery and recursive segmentation: The original sequence pair S1 and S2 are encapsulated into a root interval RootInterval, which is then delivered to the work queue as the first AnchorTask. Thread 1 executes this AnchorTask, derives the sub-interval index, searches for rare matching fragments to obtain candidate anchors, performs dynamic programming chain filtering to obtain an anchor main chain containing approximately 800 anchors, and splits the main chain to generate 801 new AnchorTasks, all of which are delivered to the work queue.

[0051] Parallel Scheduling and Execution: Threads 2 through 8 execute sub-AnchorTasks in parallel. Each thread retrieves a sub-interval task from the queue and independently performs rare anchor point search and chain filtering. If a rare anchor point is found within a sub-interval, the thread continues to split and delivers deeper AnchorTasks. If no rare anchor point is found within a sub-interval, an AlignTask is generated and delivered to the queue. The scheduling strategy prioritizes retrieving high-priority AnchorTasks from the head of the queue and allocating them to idle threads to ensure rapid recursive splitting and the generation of more parallelizable sub-tasks as quickly as possible. AlignTasks have lower priority and are only scheduled when there are no AnchorTasks to process. When a thread is idle, it retrieves an AlignTask from the queue and performs wavefront alignment. The alignment results are temporarily stored in the task node.

[0052] Results merging and performance evaluation: After the recursive task tree is formed, the main thread traverses the task tree in preorder: when it encounters an AlignTask node, it directly reads its CIGAR result; when it encounters an AnchorTask node, it reads the anchor information and recursively merges its subtask results, finally generating a complete global CIGAR string.

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

[0054] Although this paper uses terms such as input and preprocessing, rare anchor discovery, anchor chain filtering and recursive segmentation, recursive processing, and precise global comparison and completion frequently, the possibility of using other terms is not excluded. These terms are used merely to more conveniently describe and explain the essence of this invention; interpreting them as any additional limitation would contradict the spirit of this invention.

Claims

1. A rare-matching anchor-driven fast sequence alignment method for highly repetitive genomic regions, characterized in that, Includes the following steps: Step 1: Input and preprocessing: Obtain the first sequence S1 and the second sequence S2 to be compared, and perform standardization processing on S1 and S2; Step 2: Rare anchor point discovery. Construct an index structure containing the spliced ​​sequence S1 and S2. Based on the index structure, identify one or more common segments that appear at least once in both S1 and S2 and whose total number of occurrences in the spliced ​​sequence meets a preset rare threshold. The occurrence position of each common segment in S1 and S2 constitutes a candidate anchor point. Step 3: Anchor chain screening and recursive segmentation. For all candidate anchors obtained in Step 2, chain screening is performed based on collinearity constraints and a predefined scoring function to obtain a collinear anchor main chain with the highest score. Using the anchors in the anchor main chain as boundaries, S1 and S2 are segmented into multiple independent sub-interval pairs. Step 4: Recursive processing. For each sub-interval pair, repeat steps S2 and S3 until the recursion termination condition is met, resulting in multiple terminating sub-intervals. Step 5: Precise global comparison and completion. For each terminating sub-interval and the gap segment between anchor points, perform precise global comparison based on the wavefront advancement algorithm, wherein the wavefront advancement algorithm uses a segmented gap cost model for scoring. Step 6: Output the results. Combine the comparison results of all segments in Step 5 with the anchor point information in Step 3 to generate and output the complete global comparison results of S1 and S2.

2. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, The construction of the index structure and the identification of rare matching fragments in step 2 specifically include: Construct the concatenated sequence S = S1 || $ || S2, where $ is a separator that does not belong to {A,C,G,T,N}; For the spliced ​​sequence S, construct a suffix array SA, a longest common prefix array LCP, and an inverse suffix array ISA; By scanning the LCP-interval of the LCP array, common segments that appear at least once in both S1 and S2 and whose total number of occurrences in S meets the rare threshold are identified as rare matching segments.

3. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 2, characterized in that, In step 2, a rare matching fragment t is defined as a common fragment that simultaneously satisfies the following conditions: Condition 1: c_S1(t) ≥ 1 and c_S2(t) ≥ 1, where c_X(t) represents the number of times segment t appears in sequence X; Condition 2: 2 ≤ c_S(t) ≤ C_max, where c_S(t) is the total number of occurrences of segment t in the spliced ​​sequence S, and C_max is the preset maximum occurrence threshold.

4. A method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 2 or 3, characterized in that, In step 2, the search for rare matching fragments adopts a step-by-step strategy from strict to lenient: priority is given to searching for common fragments with longer lengths and fewer occurrences as candidate anchor points. When the anchor points already searched cannot cover the area to be compared, the length threshold is gradually reduced and / or the maximum occurrence threshold C_max is increased to continue the search.

5. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, The anchor chain filtering and recursive segmentation in steps 3 and 4 specifically include: Step 31: For each rare matching segment t obtained in Step 2, combine any occurrence position x of it in S1 with any occurrence position y in S2 to form an anchor point A = (x, y, w), where w = |t|, and the anchor point represents the matching correspondence between the S1 interval [x, x + w - 1] and the S2 interval [y, y + w - 1]; Step 32: Define the collinear connectable constraint: Given two anchor points Aj = (xj, yj, wj) and Ai = (xi, yi, wi), if xj + wj ≤ xi and yj + wj ≤ yi are satisfied, then Aj→Ai is said to be collinearly connectable; Step 33: Define an anchor point scoring function, including a benefit term and a penalty term: The benefit term α(i) = li / min(c_S1(ti), c_S2(ti)), where li is the length of the segment ti corresponding to the anchor point Ai, and c_X(ti) represents the number of occurrences of the segment ti in the sequence X; The penalty term β(j, i) = 0 when |(yi - yj) - (xi - xj)| = 0; otherwise β(j, i) = 2·log2(|(yi - yj) - (xi - xj)|); Step 34: Sort all candidate anchor points in ascending order of the x coordinate, and use the dynamic programming recurrence formula f(i) = max_{j < i, Aj→Ai connectable} ( f(j) + max(α(i) - β(j, i), ε) ) for chain screening to obtain the main chain of collinear anchor points with the highest score, where ε is a preset positive threshold; Step 35: Using the anchor points in the main chain of anchor points as boundaries, divide S1 and S2 into multiple mutually independent sub - interval pairs, and recursively execute Steps 2 to 35 for each sub - interval pair; Step 36: The recursive termination condition is: within the current sub - interval pair S1[l1..r1] and S2[l2..r2], no common segment that satisfies the definition of the rare matching segment in Step 2 can be found.

6. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, The segment gap cost model in Step 5 is defined as the minimum of two affine penalty functions: g_2p(k) = min( O1 + E1·(k - 1), O2 + E2·(k - 1) ), where k≥1 is the length of the continuous gap, O1, E1 are the gap opening penalty and extension penalty of the first affine penalty, O2, E2 are the gap opening penalty and extension penalty of the second affine penalty, and O2 > O1, E2 < E1.

7. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, When processing the sub - intervals generated by recursion, the local index of the sub - interval is derived by reusing the global index structure constructed in Step 2. The derivation methods include: Using the inverse suffix array ISA to extract and sort the global ranks of positions within the sub - interval to obtain the sub - interval suffix array, and calculating the sub - interval LCP array by performing range minimum query RMQ on the original LCP array.

8. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, Steps 2 to 5 are executed through a parallel computing framework, and the parallel computing framework includes: A fixed - size thread pool and a work queue; Define two task types: AnchorTask, used to execute steps 2 and 3; AlignTask, used to execute step 5. The scheduling strategy prioritizes AnchorTasks to generate more subproblems quickly. When a terminating sub-interval is generated, the corresponding AlignTask is submitted to the thread pool.

9. The method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, The standardization process in step 1 includes: converting sequence characters to uppercase, discarding non-base separators or whitespace characters, and mapping IUPAC fuzzy bases to N.

10. A method for rapid sequence alignment driven by rare matching anchors for highly repetitive genomic regions according to claim 1, characterized in that, The complete global alignment result output in step 6 includes the CIGAR string, a list of structural difference events, and the coordinate mapping relationship between S1 and S2.

Images