A haplotype-aware sequence-to-graph alignment method, system, apparatus, medium

Abstract

This invention belongs to the field of gene sequence alignment technology, specifically relating to a haplotype-aware sequence-to-graph alignment method, system, device, and medium. This method directly embeds haplotype information into the edge association structure of a pan-genome graph, enabling the graph structure itself to express haplotype paths. Then, an improved haplotype-aware partial order alignment algorithm is used to map the sequence onto this graph, allowing it to identify and utilize the embedded haplotype information, ensuring global consistency of the alignment results. Finally, a path reuse backtracking based on global caching is used to obtain the optimal alignment path, significantly improving algorithm efficiency while maintaining accuracy.

Description

Technical Field

[0001] This invention belongs to the field of gene sequence alignment technology, specifically relating to a haplotype-sensing sequence-to-image alignment method, system, device, and medium. Background Technology

[0002] High-throughput sequencing technology, with its advantages of large-scale parallelism, precise quantification, and low cost, can rapidly generate massive amounts of genomic data, becoming a core technology for analyzing the association between genetic variations and complex traits. Sequence alignment is a crucial step in high-throughput sequencing data analysis. Traditional sequence alignment compares sequencing reads with a single linear reference sequence, providing a basis for identifying variations such as single nucleotide polymorphisms and insertions / deletions. However, a single linear reference sequence cannot cover all genetic variations among individuals, especially in complex genomic regions that characterize structural variations. This limitation significantly affects the accuracy of sequence alignment in complex regions, thereby reducing the quality of variation detection.

[0003] Pangenome maps integrate the genome sequences of multiple individuals, comprehensively presenting population genetic diversity and significantly improving the detection capabilities of complex regions and structural variations. Sequence-to-map alignment, as a key technology for mapping sequencing reads to pangenome maps, utilizes multiple allelic pathways in the map to accurately locate the true source of variation in reads, which is an important prerequisite for leveraging the advantages of pangenome analysis.

[0004] However, existing sequence-to-graph alignment tools (such as partial order alignment POA) do not fully utilize haplotype path information in the graph when pursuing the optimal alignment score, ignoring global consistency at the haplotype level, resulting in deviations from the actual biological structure and insufficient biological rationality of the alignment results. At the same time, existing methods lack an effective path computation reuse mechanism, resulting in a lot of redundant computation and requiring improvement in computational efficiency. Summary of the Invention

[0005] To address the problems in the background art, the present invention provides a haplotype-sensing sequence-to-graph alignment method, system, apparatus, and medium.

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

[0007] This invention provides a haplotype-sensing sequence-to-graph alignment method, comprising:

[0008] S1: Obtain the pangenome graph, traverse all nodes in the pangenome graph, and allocate storage space for haplotypes in all predecessor edges of each node;

[0009] Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index.

[0010] S2: Based on haplotype information encoded to edge structure, perform topological sorting on the pan-genome graph, perceive the haplotype of the current path, reuse the dynamic programming state index through a caching mechanism, and minimize the total number of dynamic programming rows;

[0011] For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule.

[0012] Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score;

[0013] S3: Record the correspondence between alignment nodes and sequence positions. Before performing backtracking operations, query the global cache based on the current dynamic programming state index and use the backtracking algorithm to determine the optimal sequence to graph alignment path.

[0014] Based on the haplotype-aware sequence-to-graph alignment method described above, S2, which senses the haplotype of the current path, reuses the dynamic programming state index through a caching mechanism and minimizes the total number of dynamic programming rows, specifically:

[0015] When traversing each edge in the pangenome graph, the haplotype information stored on that edge is parsed, and a caching mechanism is maintained for each predecessor node. The caching mechanism maps the predecessor node and haplotype information to the dynamic programming state index. When the paths in the pangenome graph merge at a node, the existing dynamic programming state index is reused through the caching mechanism to generate the minimum total number of dynamic programming rows. Each dynamic programming state index is connected to the upstream dynamic programming state index of the same haplotype to obtain a consistent predecessor table.

[0016] Based on the haplotype-aware sequence-to-graph alignment method described above, the dynamic programming matrix in S2 is a state matrix with no insertion / missing penalty, a state matrix with deletion penalty, and a state matrix with insertion penalty.

[0017] Furthermore, the comparison score of the dynamic programming matrix is ​​calculated using the formula: ; ;

[0018] To achieve;

[0019] In the formula, , , These represent the alignment scores for the state matrix without insertion / missing penalty, the state matrix with deletion penalty, and the state matrix with insertion penalty, respectively. Indicates the penalty weight for haplotype inconsistency; Indicates v i The set of predecessor nodes; This represents the first step after the pangenome map is unfolded. 1 node This indicates the first sequence number of the sequencing sequence to be compared. One base, express The predecessor node, This indicates the pan-genome map node currently being compared. , This indicates the bases of the sequencing sequence currently being compared. , This indicates that the initial penalty has been deleted. This indicates a penalty for removing the extension; Indicates the initial penalty for insertion. This indicates the insertion extension penalty.

[0020] Based on the haplotype-aware sequence-to-graph alignment method described above, after calculating the alignment score of the dynamic programming matrix in step S2, the method further includes:

[0021] Record the node index corresponding to each dynamic programming state index, and associate the haplotype information of each row in the dynamic programming matrix with the corresponding node index.

[0022] Based on the haplotype-aware sequence-to-graph alignment method described above, S1 further includes storing the haplotype in bit mask form and extracting the binary bit positions through bit operations.

[0023] Based on the haplotype-aware sequence-to-graph alignment method described above, S3 operates as follows:

[0024] The CIGAR array is used to record the correspondence between alignment nodes and sequence positions, and the dynamic programming state index is converted into the corresponding node index.

[0025] Before each backtracking operation, the global cache is queried based on the current dynamic programming state index. If the path result corresponding to the dynamic programming state index exists in the global cache, the recursive backtracking calculation is skipped; if the path result corresponding to the dynamic programming state index does not exist in the global cache, the recursive backtracking algorithm is executed to generate the path result and the result is stored in the global cache.

[0026] The present invention also provides a haplotype-aware sequence-to-graph alignment system, comprising:

[0027] Haplotype encoding module: used to obtain the pan-genome graph, traverse all nodes in the pan-genome graph, and allocate storage space for the haplotypes in all predecessor edges of each node;

[0028] Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index.

[0029] Haplotype-aware partial order alignment module: Based on haplotype information encoded to edge structure, it performs topological sorting on the pan-genome graph, senses the haplotype of the current path, reuses the dynamic programming state index through a caching mechanism, and minimizes the total number of dynamic programming rows;

[0030] For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule.

[0031] Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score;

[0032] Backtracking optimization module: It is used to record the correspondence between alignment nodes and sequence positions. Before performing backtracking operation, it queries the global cache based on the current dynamic programming state index and uses the backtracking algorithm to determine the optimal sequence to graph alignment path.

[0033] The present invention also provides a haplotype-aware sequence-to-image alignment device, comprising a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the haplotype-aware sequence-to-image alignment method.

[0034] The present invention also provides a medium for storing a computer program, wherein the computer program, when executed by a processor, implements the haplotype-aware sequence-to-graph alignment method described above.

[0035] Beneficial effects

[0036] The haplotype-aware sequence-to-graph alignment method proposed in this invention directly embeds haplotype information into the edge association structure of a pan-genome graph, enabling the graph structure itself to express haplotype paths. Then, an improved haplotype-aware partial order alignment algorithm is used to map sequences onto this graph, allowing it to identify and utilize the embedded haplotype information, ensuring global consistency of the alignment results. Finally, path reuse backtracking based on global caching is used to obtain the optimal alignment path, significantly improving algorithm efficiency while maintaining accuracy. Attached Figure Description

[0037] Figure 1This is a schematic diagram of the use of haplotype information for alignment in Embodiment 1 of the present invention, wherein (a) represents sequence-to-graph alignment without haplotype sensing, and (b) represents sequence-to-graph alignment with haplotype sensing added in the present invention.

[0038] Figure 2 The diagram shows the backtracking algorithm of Embodiment 1 of the present invention, where (a) indicates no haplotype sensing and (b) indicates that haplotype sensing is added in the present invention. Detailed Implementation

[0039] The following examples are intended to illustrate the present invention, and not to further limit the invention.

[0040] Example 1

[0041] This embodiment provides a haplotype-aware sequence-to-graph alignment method, including:

[0042] S1: Obtain the pangenome graph, traverse all nodes in the pangenome graph, and allocate storage space for haplotypes in all predecessor edges of each node;

[0043] Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index.

[0044] In addition, it also includes storing the haplotype as a bitmask and extracting the binary bit positions through bit operations.

[0045] The relevant definitions in this embodiment are as follows:

[0046] A node represents a graph node in the pangenome graph, and each node stores a base sequence.

[0047] An edge represents the connection order between nodes in a pangenome graph. If there is an edge between two nodes, it means that the sequence segments represented by the two nodes are directly adjacent in at least one haplotype.

[0048] H_id represents a unique identifier for a single path, used to record which single paths pass through this edge in the graph.

[0049] dp_i represents a DP state node. During alignment, the algorithm primarily relies on dynamic programming (DP), and only physical structure nodes in the graph cannot be computed. To calculate the alignment score of the dynamic programming matrix in memory, the algorithm needs to assign an independent state index, dp_i, for each computation step. In haplotype sensing methods, the same graph node i may be visited by multiple haplotypes. In this case, the algorithm needs to expand node i into multiple independent DP states.

[0050] DP rows represent rows in the dynamic programming matrix, which are the set of score states of a base in the sequencing sequence to be aligned and all possible nodes in the graph when the algorithm processes a base.

[0051] In the specific implementation process, including single-type information storage and path calculation reuse, the specific operation of S1 is as follows:

[0052] S11: Read pangenome map in format and Sequencing sequences to be aligned in the specified format .

[0053] S12: Initialize the index structure for subsequent haplotype path management, traverse all nodes in the pangenome graph, and identify the haplotypes in all predecessor edges of each node i. Allocate storage space.

[0054] S13: Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges, and obtain the corresponding haplotype set.

[0055] S14: Stored in bitmask format. Represents 64-bit A set. Through Extract each existing The binary bit position. Incrementing count. Generate a globally unique dynamic programming state index Each Allocated the only Establish predecessor nodes arrive Then The mapping.

[0056] The above haplotype information storage process is implemented by the algorithm in Table 1.

[0057] Table 1. Haplotype Information Storage Process

[0058]

[0059] Within each edge, the haplotype information storage algorithm iterates through all haplotypes that pass through that edge. The sum of this operation equals the sum of the number of times all haplotypes are traversed along that edge. The time complexity is O(n log n). Where H represents the number of haplotypes. L nodes This represents the number of times a haplotype is recorded during graph traversal. There are H distinct haplotypes in the graph, and traversing the graph requires... L nodes If there are 1 edge, then the time complexity of the single-type information storage algorithm is O(H×). L nodes ).

[0060] S15: Add a new cache structure to record... arrive The mapping. During subsequent comparisons, the same mapping was found in the cache. It can be reused directly without secondary calculation.

[0061] The above path calculation reuse process is implemented by the algorithm in Table 2.

[0062] Table 2 Path Calculation Reuse Process

[0063]

[0064] A new dynamic programming line (dp_i) is created only when a haplotype first appears on a node (node_id); otherwise, all subsequent paths encountering the (node_id, H_id) combination will reuse the already created DP line, i.e., the dynamic programming line.

[0065] Through the above operations, this invention embeds haplotype information into the pan-genome graph structure and establishes a path computation reuse mechanism to reduce redundant computation and improve algorithm efficiency.

[0066] S2: Based on haplotype information encoded to edge structure, perform topological sorting on the pan-genome graph, perceive the haplotype of the current path, reuse the dynamic programming state index through a caching mechanism, and minimize the total number of dynamic programming rows;

[0067] For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule.

[0068] Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score.

[0069] In the specific implementation process, including dynamic programming based on haplotype perception, S2 operates as follows:

[0070] S21: Topological sorting of the pangenome graph based on haplotype information encoded to edge structure.

[0071] S22: Aware of the singlet nature of the current path, reuse the dynamic programming state index through a caching mechanism and minimize the total number of dynamic programming rows, specifically:

[0072] When traversing each edge in the pangenome graph, the haplotype information stored on that edge is parsed for each predecessor node. Maintain a caching mechanism to link predecessor nodes with haplotype information. Mapping to dynamic programming state index When paths in the pangenome graph converge at nodes, existing dynamic programming state indexes are reused through a caching mechanism. Generate the minimum total number of rows in the dynamic programming. Index each dynamic programming state Upstream dynamic programming state index with the same haplotype Connect them to obtain a consistent predecessor table.

[0073] S23: For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule.

[0074] The dynamic programming matrix is ​​a state matrix with no insertion / missing penalty, a state matrix with deletion penalty, and a state matrix with insertion penalty.

[0075] Preferably, the comparison score of the dynamic programming matrix is ​​calculated using the formula: ; ;

[0076] To achieve this.

[0077] In the formula, , , These represent the alignment scores for the state matrix without insertion / missing penalty, the state matrix with deletion penalty, and the state matrix with insertion penalty, respectively. Indicates the penalty weight for haplotype inconsistency; Indicates v i The set of predecessor nodes; This represents the first step after the pangenome map is unfolded. 1 node This indicates the first sequence number of the sequencing sequence to be compared. One base, express The predecessor node, This indicates the pan-genome map node currently being compared. , This indicates the bases of the sequencing sequence currently being compared. , This indicates that the initial penalty has been deleted. This indicates a penalty for removing the extension; Indicates the initial penalty for insertion. This indicates the insertion extension penalty.

[0078] In sequence-to-graph alignment, the insertion penalty state matrix represents the state of graph node v when an insertion is introduced into sequence q. i Since it remains unchanged, the Y-state matrix does not involve predecessor nodes in the graph structure, that is, there is no p∈pred(v i The insertion operation corresponds to skipping the current base q in the sequence during the alignment process. j Therefore, its state transition should come from the previous column, and the graph nodes remain unchanged, i.e. (i,j-1).

[0079] The dynamic programming process based on haplotype sensing described above is implemented by the algorithm in Table 3, with a time complexity of O(n log n). Where P represents the total number of rows in the DP matrix, and Q represents the length of the query sequence, i.e., the total number of bases in the sequencing sequence to be aligned. In haplotype-aware sequence-to-graph alignment, during dynamic programming, a graph node is expanded into multiple independent one-to-one DP states based on the multiple haplotypes that pass through it. P is the total number of all expanded DP states in the graph.

[0080] To find the optimal scoring path, the algorithm traverses each DP state along the graph's topology, performing a total of P loops. Within each loop, the algorithm calculates the alignment score between the current graph state and the entire query sequence. The final time complexity is P states multiplied by Q bases.

[0081] Table 3. Dynamic Programming Process Based on Haplotype Sensing

[0082]

[0083] S24: After calculating the comparison score of the dynamic programming matrix, the following steps are also included:

[0084] pass Record the node index i corresponding to each dynamic programming state index, and associate the haplotype information of each row in the dynamic programming matrix with the corresponding node index.

[0085] S25: To obtain the optimal haplotype sensing alignment path, traverse all predecessor edges of each node and find the corresponding nodes in the current graph. Corresponding haplotype By matching the same Determine the haplotype continuity path.

[0086] S26: Calculate the best alignment score. The procedure is as follows:

[0087] The columns of the dynamic programming matrix represent the base positions of the query sequence. To align the entire query sequence onto the graph, the optimal alignment score must be in the last column of the matrix, which is the position with coordinate Q, indicating that the entire query sequence has been aligned.

[0088] The comparison must end at the end of the query sequence, not necessarily at the last node in the graph; a virtual convergence node needs to be set. All graph nodes that are allowed as comparison endpoints will have a logical edge connected to the convergence node.

[0089] The algorithm queries the last column of the dynamic programming matrix to find the comparison scores of all legal exit nodes. By traversing these final scores, the maximum value is selected as the optimal comparison score.

[0090] When advancing the dynamic programming matrix, if the path crosses different haplotypes, a haplotype inconsistency penalty weight s is deducted from the score. If the path follows the same haplotype, there is no additional penalty. This step achieves haplotype awareness.

[0091] Through the above operations, this invention improves the traditional partial-order alignment algorithm. It can perceive and utilize haplotype information during the dynamic programming process. By introducing a haplotype inconsistency penalty, it calculates the optimal alignment score that is more in line with biological significance, thus ensuring the biological rationality of the alignment results.

[0092] besides, Figure 1 The diagram illustrates the process of sequence-to-graph alignment without haplotype sensing and the alignment process of incorporating the haplotype sensing algorithm of this invention. The three haplotypes are represented by different colors, and the alignment path is indicated by a red dashed line.

[0093] Both studies used the same pan-genome directed acyclic graph and query sequence as input, but obtained different alignment results.

[0094] Figure 1In (a), it can be observed that sequence-to-graph alignment without haplotype awareness typically pursues the highest base match score, regardless of whether the alignment path actually exists, as long as the path exists in the graph. The final result is a chimeric path pieced together from variations of different individuals. Although it achieves the optimal score, it is highly likely an illegal path that does not exist biologically.

[0095] Figure 1 In (b), haplotype-aware sequence-to-graph alignment was incorporated. The algorithm no longer simply pursues the highest score at the base level, but introduces a haplotype continuity constraint. If an alignment path attempts to jump from one haplotype to another, a penalty point is deducted. The algorithm ultimately outputs a coherent haplotype path with real biological significance, avoiding the generation of chimeric paths.

[0096] This invention strikes a balance between alignment scores and biological continuity, ensuring that alignment results are biologically plausible and interpretable while maintaining alignment accuracy. This is crucial for accurate variant recall and biological function analysis using pangenomes.

[0097] S3: Record the correspondence between alignment nodes and sequence positions. Before performing backtracking, query the global cache based on the current dynamic programming state index and use the backtracking algorithm to determine the optimal sequence-to-graph alignment path. The specific operations are as follows:

[0098] The CIGAR array is used to record the correspondence between alignment nodes and sequence positions, and the dynamic programming state index is converted into the corresponding node index.

[0099] Before each backtracking operation, the global cache is queried based on the current dynamic programming state index. If the path result corresponding to the dynamic programming state index exists in the global cache, the recursive backtracking calculation is skipped; if the path result corresponding to the dynamic programming state index does not exist in the global cache, the recursive backtracking algorithm is executed to generate the path result and the result is stored in the global cache.

[0100] The above-mentioned path reuse backtracking process based on global caching is implemented by the algorithm in Table 4.

[0101] Table 4. Path Reuse Backtracking Algorithm Based on Global Caching

[0102]

[0103] This invention relates to a path reuse backtracking algorithm based on global caching, which relies on the management of haplotype information. The caching ensures the uniqueness of path segments and the safety of reuse.

[0104] Figure 2The diagram illustrates the construction and backtracking process of the score matrix for sequence-to-graph alignment without haplotype sensing compared to the sequence-to-graph alignment with haplotype sensing incorporated in this invention. The query sequence serves as the rows of the matrix, and the pan-genome graph as the columns. It can be seen that, compared to... Figure 2 Compared to the sequence-to-graph alignment process without haplotype sensing in (a), after adding haplotype sensing, Figure 2 In (b), the backtracking was strictly performed according to the haplotype path, which ensured the consistency of biology.

[0105] To further verify the effectiveness of the present invention, the following experiments were conducted:

[0106] A 5000bp reference sequence was randomly selected from the human reference genome GRCh38. (Using...) Ten simulated sequences were generated for each structural mutation rate of 0.1%, 1%, 1.5%, and 2%. Using simulated sequences with different mutation rates aimed to test whether the algorithm could maintain high alignment capability when faced with increasingly complex graph structures.

[0107] Subsequently, a pan-genome map was constructed by combining these reference sequences and simulated sequences, and then... After generating a simulated dataset with a depth of 1000, negative links were filtered to test the accuracy of HaPOA.

[0108] HaPOA is a sequence-to-graph alignment tool with haplotype sensing designed based on the haplotype sensing algorithm and the sequence-to-graph alignment tool abPOA.

[0109] KSW stands for Striped Smith-Waltman Algorithm, specifically designed for calculating dynamic programming score matrices. It natively supports complex affine vacancy penalties, making it suitable for handling insertions, missing values, and mutations. KSW utilizes Single Instruction Multiple Data Stream (SIMD) acceleration, allowing it to simultaneously calculate scores for multiple cells in a dynamic programming problem.

[0110] The main differences between the comparison method ksw and our method HaPOA are:

[0111] In the dynamic programming state space, one physical graph node in this method corresponds to multiple haplotype DP states, while in the comparative method, one base corresponds to one DP state.

[0112] In the scoring mechanism, this method incorporates a penalty for inconsistent haplotype scores, which is not present in the comparison method.

[0113] In the backtracking process, the comparison method only backtracks based on the matrix direction with the higher score. In this method, the backtracking is based on the predecessor edge and the haplotype consistency check.

[0114] In a biological sense, the contrastive method only looks for the highest-scoring path, and the path found is likely to be a chimeric path, that is, a path that does not exist in itself. However, this method incorporates haplotype sensing to find real continuous sequences, which is more in line with biological consistency.

[0115] Table 5 Comparison of memory usage and runtime between ksw2 and HaPOA

[0116]

[0117] According to the results in Table 5, HaPOA demonstrated a significant advantage in genome alignment tasks, with a runtime of 952.588 seconds, a 40.3% improvement over ksw's 1595.716 seconds. This indicates that HaPOA significantly improved alignment speed through graph structure, single-instruction multiple-data stream (SIMD) acceleration, and path reuse. However, HaPOA requires traversing every edge in the graph and completely saving haplotype paths, resulting in a peak memory usage of 18.294 GB. Comparative experiments further validated that both HaPOA and ksw achieved 100% accuracy in finding correct haplotype information during testing, proving that HaPOA maintains both high speed and accuracy in correctly aligned haplotypes.

[0118] In summary, this invention proposes a haplotype information storage algorithm, a haplotype-aware dynamic programming algorithm, and a haplotype-aware path reuse calculation backtracking algorithm by embedding haplotype information into the edge structure, thereby realizing a fast and accurate method for haplotype-aware sequence to graph alignment.

[0119] The haplotype-aware sequence-to-graph alignment method proposed in this invention directly embeds haplotype information into the edge association structure of a pan-genome graph, enabling the graph structure itself to express haplotype paths. Then, an improved haplotype-aware partial order alignment algorithm is used to map sequences onto this graph, allowing it to identify and utilize the embedded haplotype information, ensuring global consistency of the alignment results. Finally, path reuse backtracking based on global caching is used to obtain the optimal alignment path, significantly improving algorithm efficiency while maintaining accuracy.

[0120] The present invention also provides a haplotype-aware sequence-to-graph alignment system, comprising:

[0121] Haplotype encoding module: used to obtain the pan-genome graph, traverse all nodes in the pan-genome graph, and allocate storage space for the haplotypes in all predecessor edges of each node;

[0122] Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index.

[0123] Haplotype-aware partial order alignment module: Based on haplotype information encoded to edge structure, it performs topological sorting on the pan-genome graph, senses the haplotype of the current path, reuses the dynamic programming state index through a caching mechanism, and minimizes the total number of dynamic programming rows;

[0124] For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule.

[0125] Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score;

[0126] Backtracking optimization module: It is used to record the correspondence between alignment nodes and sequence positions. Before performing backtracking operation, it queries the global cache based on the current dynamic programming state index and uses the backtracking algorithm to determine the optimal sequence to graph alignment path.

[0127] The present invention also provides a haplotype-aware sequence-to-image alignment device, comprising a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the haplotype-aware sequence-to-image alignment method.

[0128] The present invention also provides a medium for storing a computer program, wherein the computer program, when executed by a processor, implements the haplotype-aware sequence-to-graph alignment method described above.

Claims

1. A haplotype-sensing sequence-to-graph alignment method, characterized in that, include: S1: Obtain the pangenome graph, traverse all nodes in the pangenome graph, and allocate storage space for haplotypes in all predecessor edges of each node; Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index. S2: Based on haplotype information encoded to edge structure, perform topological sorting on the pan-genome graph, perceive the haplotype of the current path, reuse the dynamic programming state index through a caching mechanism, and minimize the total number of dynamic programming rows; For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule. Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score; S3: Record the correspondence between alignment nodes and sequence positions. Before performing backtracking operations, query the global cache based on the current dynamic programming state index and use the backtracking algorithm to determine the optimal sequence to graph alignment path.

2. The haplotype-sensing sequence-to-graph alignment method according to claim 1, characterized in that, S2 describes the perception of the current path's singlet form, reusing the dynamic programming state index through a caching mechanism and minimizing the total number of dynamic programming rows, specifically as follows: When traversing each edge in the pangenome graph, the haplotype information stored on that edge is parsed, and a caching mechanism is maintained for each predecessor node. The caching mechanism maps the predecessor node and haplotype information to the dynamic programming state index. When the paths in the pangenome graph merge at a node, the existing dynamic programming state index is reused through the caching mechanism to generate the minimum total number of dynamic programming rows. Each dynamic programming state index is connected to the upstream dynamic programming state index of the same haplotype to obtain a consistent predecessor table.

3. The haplotype-sensing sequence-to-graph alignment method according to claim 1, characterized in that, The dynamic programming matrix described in S2 is a state matrix with no insertion / missing penalty, a state matrix with deletion penalty, and a state matrix with insertion penalty.

4. The haplotype-sensing sequence-to-graph alignment method according to claim 3, characterized in that, The comparison score of the dynamic programming matrix is ​​calculated using the formula: ； ； To achieve; In the formula, , , These represent the alignment scores for the state matrix without insertion / missing penalty, the state matrix with deletion penalty, and the state matrix with insertion penalty, respectively. Indicates the penalty weight for haplotype inconsistency; Indicates v i The set of predecessor nodes; This represents the first step after the pangenome map is unfolded. One node; This indicates the first sequence number of the sequencing sequence to be compared. One base; express The predecessor node; This indicates the pangenome map node currently being compared. , This indicates the bases of the sequencing sequence currently being compared. , This indicates that the initial penalty has been deleted. This indicates a penalty for removing the extension; Indicates the initial penalty for insertion. This indicates the insertion extension penalty.

5. The haplotype-sensing sequence-to-graph alignment method according to claim 1, characterized in that, After calculating the comparison score of the dynamic programming matrix as described in S2, the method further includes: Record the node index corresponding to each dynamic programming state index, and associate the haplotype information of each row in the dynamic programming matrix with the corresponding node index.

6. The haplotype-sensing sequence-to-graph alignment method according to claim 1, characterized in that, S1 also includes storing the haplotype as a bit mask and extracting the binary bit positions through bit operations.

7. The haplotype-sensing sequence-to-graph alignment method according to claim 1, characterized in that, The specific operation of S3 is as follows: The CIGAR array is used to record the correspondence between alignment nodes and sequence positions, and the dynamic programming state index is converted into the corresponding node index. Before each backtracking operation, the global cache is queried based on the current dynamic programming state index. If the path result corresponding to the dynamic programming state index exists in the global cache, the recursive backtracking calculation is skipped; if the path result corresponding to the dynamic programming state index does not exist in the global cache, the recursive backtracking algorithm is executed to generate the path result and the result is stored in the global cache.

8. A haplotype-sensing sequence-to-graph alignment system, characterized in that, include: Haplotype encoding module: used to obtain the pan-genome graph, traverse all nodes in the pan-genome graph, and allocate storage space for the haplotypes in all predecessor edges of each node; Traverse all outgoing edges of the predecessor node of each node, filter out the outgoing edges pointing to the current node, extract the haplotype identifier carried in the outgoing edges to obtain the corresponding haplotype set, assign a unique dynamic programming state index to each haplotype, and establish a mapping between predecessor node-haplotype-dynamic programming state index. Haplotype-aware partial order alignment module: Based on haplotype information encoded to edge structure, it performs topological sorting on the pan-genome graph, senses the haplotype of the current path, reuses the dynamic programming state index through a caching mechanism, and minimizes the total number of dynamic programming rows; For the sequencing sequences to be aligned, the alignment score of the dynamic programming matrix is ​​calculated by introducing haplotype inconsistency penalty weights according to the dynamic programming state index order and the matrix penalty rule. Traverse all predecessor edges of each node, find the haplotype corresponding to the current node, determine the haplotype continuity path, and calculate the best alignment score; Backtracking optimization module: It is used to record the correspondence between alignment nodes and sequence positions. Before performing backtracking operation, it queries the global cache based on the current dynamic programming state index and uses the backtracking algorithm to determine the optimal sequence to graph alignment path.

9. A haplotype-sensing sequence-to-image alignment device, characterized in that, It includes a processor and a memory, wherein the processor executes a computer program stored in the memory to implement a haplotype-aware sequence-to-graph alignment method as described in any one of claims 1-7.

10. A medium, characterized in that, Used to store a computer program, wherein the computer program, when executed by a processor, implements a haplotype-aware sequence-to-graph alignment method as described in any one of claims 1-7.

Images