A shingled spectral hidden Markov sequence alignment method on a directed acyclic pan-genome graph

Abstract

The invention discloses a shingled spectral hidden Markov sequence alignment method on a directed acyclic pan-genome graph in the technical field of sequence alignment. The method comprises the following steps: pre-processing the directed acyclic pan-genome graph before alignment to obtain a Viterbi graph, a training graph and an alignment reference graph, then using the length of the longest alignment reference path in the alignment reference graph as the number of matching states of an shingled spectral hidden Markov model, training on the training graph after initializing the model parameters and the shingled width of each node to obtain an shingled spectral hidden Markov model, and when multiple sequence alignment is required, using a Viterbi decoding algorithm based on virtual node probability calculation on the model to calculate the most likely state path of each sequence in the shingled spectral hidden Markov model on the Viterbi graph, determining the alignment relationship of each position in the sequence according to the state path, thereby completing the multiple sequence alignment. The invention greatly reduces the calculation time and storage space requirements of the sequence alignment.

Description

Technical Field

[0001] The invention relates to the technical field of sequence alignment, in particular to a shingled hidden Markov sequence alignment method on a directed acyclic pan-genome graph. Background Art

[0002] With the rapid development of sequencing technology in recent years, massive genomic data sets have emerged. Since these genomes contain a large amount of information that records the evolutionary process and is crucial for our understanding and manipulation of synthetic biological systems, reliable analysis and comparison of such data sets naturally become a rigid demand. For example, the number of sequences of the new coronavirus genome has exceeded 20 million, and the number of sequences of the human genome is expected to exceed 100 million in the next few years. Many other important species (pathogenic microorganisms, pets, crops, etc.) have also begun to be sequenced in large quantities. One of the most important foundations of biological sequence analysis is multiple sequence alignment.

[0003] There are many different multiple sequence alignment techniques. They can be roughly divided into three categories: 1) various alignment techniques based on dynamic programming algorithms; 2) alignment techniques based on spectral hidden Markov models; 3) hybrid techniques of these two techniques. These techniques have been around for many years and have made great contributions to people's understanding of biological sequences.

[0004] However, none of the above sequence alignment technologies can effectively cope with high-precision multi-sequence alignment of massive genomes. The reason is that the computing power required for these technologies to perform high-precision alignment of sequence data sets of more than 100,000 is far beyond the budget of general research groups and companies. The computing power and computing time required for sequence data sets of more than 10 million are beyond the reach of large research institutions and companies. Taking 4 million high-quality new coronavirus genomes as an example, if consistency-based high-precision alignment technology is used, it will take more than 1,000 years of calculation on the fastest E-class supercomputer currently available. Summary of the invention

[0005] The purpose of this section is to summarize some aspects of the embodiments of the present invention and briefly introduce some preferred embodiments. Some simplifications or omissions may be made in this section and the specification abstract and the invention title of this application to avoid blurring the purpose of this section, the specification abstract and the invention title, and such simplifications or omissions cannot be used to limit the scope of the present invention.

[0006] Therefore, the purpose of the present invention is to provide a shingled hidden Markov sequence alignment method on a directed acyclic pan-genome graph, which greatly reduces the computation time and storage space requirements.

[0007] To solve the above technical problems, according to one aspect of the present invention, the present invention provides the following technical solutions:

[0008] The steps of the shingled hidden Markov sequence alignment method on the directed acyclic pan-genome graph are as follows:

[0009] S1. The directed acyclic pan-genome graph is subjected to pre-alignment pre-processing by reducing the fragment length, removing degenerate base nodes and low-weight nodes, and obtaining a Viterbi graph, a training graph and an alignment reference graph;

[0010] S2, taking the length of the longest comparison reference path in the comparison reference graph as the number of matching states of the shingled spectral hidden Markov model, and performing training on the training graph after initializing the model parameters and the shingled width of each node to obtain the shingled spectral hidden Markov model;

[0011] S3. After obtaining the shingled spectral hidden Markov model, the Viterbi decoding algorithm based on virtual node probability calculation is used to calculate the most likely state path of each sequence in the shingled spectral hidden Markov model on the Viterbi graph, and the alignment relationship of each position in the sequence is determined according to the state path, thereby completing the multiple sequence alignment.

[0012] As a preferred solution of the shingled hidden Markov sequence alignment method on the directed acyclic pan-genome graph of the present invention, in step S1, the step of reducing the length of the directed acyclic pan-genome graph segment is as follows:

[0013] Decompose all starting nodes of the directed acyclic pan-genome graph into L1-L2+1 new nodes, and each node retains a sequence segment of length L2;

[0014] For all other nodes, only the last L2 characters are retained in the sequence fragments;

[0015] Node merging based on topological sequence coordinates is performed in the graph, wherein no node merging is performed when the length of the last segment is reduced to 1 in the stepwise reduction of the segment length.

[0016] As a preferred scheme of the shingled spectral hidden Markov sequence alignment method on the directed acyclic pan-genome graph described in the present invention, in step S1, the step of degenerate base nodes of the directed acyclic pan-genome graph is: all nodes and corresponding edges containing degenerate bases in the directed acyclic pan-genome graph are removed.

[0017] As a preferred solution of the shingled hidden Markov sequence alignment method on the directed acyclic pan-genome graph of the present invention, in step S1, the step of removing low-weight nodes of the directed acyclic pan-genome graph is as follows:

[0018] First, all nodes below the global low weight threshold are removed from the directed acyclic pan-genome graph;

[0019] Then, traverse from each tail node in reverse topology order, stop traversing at the first node whose weight is greater than the tail low weight threshold, and delete the tail low weight nodes that have been traversed;

[0020] Among them, the global low weight threshold is set to 0.01 or 0.001, and the tail low weight threshold is set to 0.01.

[0021] As a preferred embodiment of the shingled hidden Markov sequence alignment method on the directed acyclic pan-genome graph of the present invention, in step S2, the shingled hidden Markov model uses three different smoothing parameters of exp(-3), exp(-5) and exp(-7), and these specific values ​​are for reference. For simplicity, the same value can be used for all probability variables, and obviously different values ​​can also be used. Since the high-dimensional space optimization problem is closely related to the starting point, people can use various random sampling algorithms according to actual conditions to further optimize the initialization value of the parameter, and the formula used for smoothing is as follows:

[0022]

[0023] Among them, for a selected set of n possible variables, p i is the probability of the ith possibility, p j is the probability of the jth possibility, c i is the smoothing parameter for the ith possibility, c j is the smoothing parameter for the jth possibility, i is the updated probability.

[0024] As a preferred embodiment of the shingled hidden Markov model alignment method on the directed acyclic pan-genome graph of the present invention, in step S2, the parameters of the shingled hidden Markov model include the starting probability of the starting state, the transition probabilities from the starting state to insertion (I), deletion (D) and matching (M) are 1 / 3 respectively, the emission probabilities of the insertion parts are initialized to 1 / 4, and the transition probabilities of all matching state positions are given the same initialization value, wherein T MM = 0.96, T MI = T MD = exp(-4), T DI = T ID = 0, T II = T IM = T DM = T DD = 0.5, T M-end= 0.5. Of course, these specific values ​​are for reference only. For simplicity, the same value can be used for all matching state positions. Obviously, different values ​​can also be used. Since the high-dimensional space optimization problem is closely related to the starting point, people can use various random sampling algorithms to further optimize the initialization values ​​of the parameters according to the actual situation. In addition, the emission probability of each matching state is initialized to the proportion of various bases in the corresponding topological sequence coordinate node.

[0025] As a preferred embodiment of the shingled hidden Markov model sequence alignment method on the directed acyclic pan-genome graph of the present invention, in step S2, the shingled hidden Markov model performs the following steps on the input directed acyclic pan-genome graph with alignment:

[0026] First, we traverse the order-complement dependency, treat each path without branches in the directed acyclic pan-genome graph as a separate computing task, and record the dependencies between different computing tasks. All computing tasks that meet the conditions are placed in a task queue.

[0027] The idle thread extracts the task from the task queue, and after completion, it records and checks whether the subsequent tasks that depend on it can be executed. If all dependencies are met, they are added to the tail of the task queue.

[0028] As a preferred scheme of the shingled spectral hidden Markov sequence alignment method on the directed acyclic pan-genome graph described in the present invention, in the shingled spectral hidden Markov model, a single forward virtual parent / child node is constructed for the graph node with multiple parent / child nodes in the training data set and the probability summation calculation is performed, and a single backward virtual parent / child node is constructed and the probability summation calculation is performed.

[0029] As a preferred solution of the shingled spectral hidden Markov sequence alignment method on the directed acyclic pan-genome graph described in the present invention, in order to realize the parallel training of the shingled spectral hidden Markov model, each branchless path on the directed acyclic graph is taken as a computing task, and the dependency between these tasks depends on the topological order between them. Tasks that are not dependent on each other are calculated simultaneously using different threads.

[0030] As a preferred scheme of the shingled spectral hidden Markov sequence alignment method on the directed acyclic pan-genome graph described in the present invention, in the Viterbi decoding algorithm of the shingled spectral hidden Markov model, a virtual parent node i-1 is constructed for the graph node with multiple parent / child nodes in the training data set to implement the corresponding forward maximum probability calculation, and in the backtracking stage of the Viterbi decoding algorithm, in the case of multiple parent nodes, only one parent node path is determined as the maximum probability path in the initial backtracking, and the paths corresponding to the other parent nodes are backtracked from the child nodes again. In the case of multiple child nodes, the virtual path is verified and excluded, and the possible different states of different genomes in the same node are recorded.

[0031] Compared with the prior art, the present invention has the following beneficial effects:

[0032] 1. In the shingled spectral hidden Markov model of the present invention, a single virtual parent / child node is constructed for a graph node with multiple parent / child nodes in a directed acyclic pan-genome graph, thereby smoothly realizing forward and backward calculations in hidden Markov model training, achieving a high-speed calculation rate, and improving the comparison rate.

[0033] 2. The directed acyclic pan-genome graph of the present invention is subjected to pre-alignment preprocessing by reducing the fragment length, removing degenerate base nodes and low-weight nodes, which greatly reduces the complexity of the pan-genome graph while maintaining the effective information in the graph, and improves the processing accuracy and speed of the shingled spectral hidden Markov model.

[0034] 3. The present invention uses a shingled hidden Markov model to traverse the input directed acyclic pan-genome graph with alignment, and treats each path without branches in the graph as a separate computing task, and records the dependencies between different computing tasks. All computing tasks that meet the conditions are placed in a task queue. The idle thread extracts tasks from the task queue, and after completion, records and checks whether the subsequent tasks that depend on it can be executed. If all dependencies are met, they are added to the tail of the task queue for processing, so as to realize parallel computing and further improve the processing speed.

[0035] 4. The present invention aims at the Viterbi decoding algorithm, and constructs a virtual parent node i-1 for the graph node with multiple parent / child nodes in the training data set to implement the corresponding forward maximum probability calculation.

[0036] 5. The present invention adopts shingled training for the shingled spectral hidden Markov model, which greatly reduces the calculation time and storage space requirements compared with the ordinary spectral hidden Markov model. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the present invention will be described in detail below in combination with the accompanying drawings and detailed embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative labor. Among them:

[0038] Figure 1 A directed acyclic pan-genome graph alignment preprocessing flow chart provided by the present invention;

[0039] Figure 2 A schematic diagram of converting a linear path in a directed acyclic pan-genome graph provided by the present invention into a parallel computing task;

[0040] Figure 3 Schematic diagram of the construction of the virtual parent node for forward calculation and its probability summation provided by the present invention, wherein (a) is a virtual parent node with m i Node i of a parent node, (b) is a schematic diagram of the virtual parent node (i-1), (c) the sum of the probabilities of my virtual parent node, (d) is a schematic diagram of the calculation from the virtual parent node to node i, and (e) is the formula for obtaining the probability of a virtual parent node by summing the probabilities of multiple parent nodes;

[0041] Figure 4 A schematic diagram of the construction of the backward computing virtual subnodes and their probability summation provided by the present invention, wherein (a) there are n i Node i of 1 child nodes, (b) is a schematic diagram of the virtual child node (i+1), (c) is the sum of the probabilities of the virtual child nodes, and (d) is a schematic diagram of the calculation of the probability from the virtual child node to node i;

[0042] Figure 5 Schematic diagram of the correspondence between the shingled hidden Markov model and the comparison reference graph and the reference path topology sequence coordinate range of the non-comparison reference path nodes, where (a) is a schematic diagram of the correspondence between the comparison reference graph (bottom) and the shingled hidden Markov model (top), and the number of matching states is consistent with the length of the comparison reference path. The nodes on the comparison reference path are displayed in orange, and the non-reference path nodes are represented by gray dotted circles. The number in the square above each reference path node represents the weight of the node, (b) non-reference path nodes are taken into consideration, represented by solid circles, (c) the reference path topology sequence coordinate range corresponding to the green "T" node is "1-4", (d) the reference path topology sequence coordinate range corresponding to the green "C" node is "7-10", (e) the reference path topology sequence coordinate range corresponding to the green "T" node at the bottom of the figure is "5-9", (f) is a schematic diagram of the reference path topology sequence coordinate range of all non-reference path nodes being established;

[0043] Figure 6The present invention provides a schematic diagram of constructing a virtual parent node (i-1) to realize the corresponding forward maximum probability calculation, wherein a) is m i parent nodes, b) is a schematic diagram of the virtual child node (i-1), c) is the sum of the probabilities of the virtual parent nodes, and d) is a schematic diagram of the calculation of the probability from the virtual parent node (i-1) to the node i. DETAILED DESCRIPTION

[0044] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0045] Secondly, the present invention is described in detail with reference to schematic diagrams. When describing the embodiments of the present invention in detail, for the sake of convenience, the cross-sectional diagrams showing the device structure will not be partially enlarged according to the general scale, and the schematic diagrams are only examples, which should not limit the scope of protection of the present invention. In addition, in actual production, the three-dimensional dimensions of length, width and depth should be included.

[0046] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0047] The present invention provides a shingled hidden Markov sequence alignment method on a directed acyclic pan-genome graph, and the specific steps are as follows:

[0048] S1. The directed acyclic pan-genome graph is subjected to pre-alignment pre-processing by reducing the fragment length, removing degenerate base nodes and low-weight nodes, and obtaining a Viterbi graph, a training graph and an alignment reference graph;

[0049] S2, taking the length of the longest comparison reference path in the comparison reference graph as the number of matching states of the shingled spectral hidden Markov model, and performing training on the training graph after initializing the model parameters and the shingled width of each node to obtain the shingled spectral hidden Markov model;

[0050] S3. After obtaining the shingled spectral hidden Markov model, the Viterbi decoding algorithm based on virtual node probability calculation is used to calculate the most likely state path of each sequence in the shingled spectral hidden Markov model on the Viterbi graph, and the alignment relationship of each position in the sequence is determined according to the state path, thereby completing the multiple sequence alignment.

[0051] The pre-alignment pre-processing of the directed acyclic pan-genome graph in step S1 above, which includes reducing the fragment length, removing degenerate base nodes and low-weight nodes, is specifically explained as follows:

[0052] At present, the construction technology based on topological sequence is the only efficient and feasible method for constructing a directed acyclic pan-genome graph. In order to avoid too many topological sequence conflict calculations during the construction process, the construction of this graph is carried out in units of longer sequence fragments. However, in spectral hidden Markov alignment, it is more convenient to represent the graph in units of single bases / amino acids (or other related characters). To achieve this goal, we have developed a step-by-step reduction algorithm for fragment length, so that the fragment length in the directed acyclic graph based on fragment topological sequence (FTO-DAG) can be reduced to 1, and the graph representation using single characters can be realized. In addition, some nodes / paths in the pan-genome graph will contain degenerate characters (that is, characters that cannot be determined with a high probability, and are represented by different alternative characters according to the possible character value range, such as "N" represents any character in "A, T, C, G"). On the one hand, the nodes containing these bases increase the complexity and training cost of the pan-genome graph, but on the other hand, they cannot provide more effective information. Therefore, these nodes can be removed. When the shingled hidden Markov model is used to process the directed acyclic pan-genome graph, it does not significantly affect the model accuracy. Finally, some low-weight nodes in the graph may be sequencing errors, so removing these low-weight nodes can significantly reduce the complexity of the pan-genome graph while maintaining the effective information in the graph. The overall pre-processing scheme of the directed acyclic pan-genome graph before alignment is as follows: Figure 1 shown.

[0053] 1. FTO-DAG fragment length reduction technology: that is, to change the FTO-DAG with a fragment length of L1 into a FTO-DAG with a fragment length of L2 (L1>L2), the steps are as follows:

[0054] 1) Decompose all starting nodes into (L1-L2+1) new nodes, and each node retains a sequence segment of length L2;

[0055] 2) Sequence fragments in all other nodes only retain the last L2 characters;

[0056] 3) Nodes are merged in the graph based on topological sequence coordinates. However, in order to retain sufficient adjustment space for the alignment of character positions of different sequences, no node merging is performed when the length of the last segment is reduced to one (a process we call the atomization of the directed acyclic pan-genome graph), and the penultimate segment length determines the flexibility of the adjustment of the character positions of different sequences in the alignment.

[0057] 2. Degenerate base node removal: All nodes containing degenerate bases and the corresponding edges are removed from the directed acyclic pan-genome graph.

[0058] 3. Removal of low-weight nodes: Different locations of the genome have different degrees of conservation and sequencing errors. Therefore, different standards can be used to define low-weight nodes for different regions. For the new coronavirus genome, we used different low-weight definition thresholds for the tail of the genome and other regions, namely the global low-weight threshold and the tail low-weight threshold. First, remove all nodes below the global low-weight threshold from the directed acyclic pan-genome graph, and then start traversing from each tail node in reverse topology order. Stop traversing at the first node whose weight is greater than the tail low-weight threshold, and delete the tail low-weight nodes that have been traversed.

[0059] For step S2, the shingled spectrum hidden Markov model is explained in detail as follows:

[0060] 1. The shingled spectral hidden Markov model uses three different smoothing parameters, exp(-3), exp(-5) and exp(-7), to eliminate the possibility that some possible options are completely abandoned due to underflow due to low probability during the calculation process. Specifically, first remove all nodes below the global low weight threshold from the directed acyclic pan-genome graph, then start traversing from each tail node in reverse topological order, stop traversing at the first node whose weight is greater than the tail low weight threshold, and delete the tail low weight nodes that have been traversed. The smoothing formula used is as follows:

[0061]

[0062] Among them, for a selected set of n possible variables, p i is the probability of the ith possibility, p j is the probability of the jth possibility, c i is the smoothing parameter for the ith possibility, c j is the smoothing parameter for the jth possibility, i is the updated probability.

[0063] The parameters of the shingled spectral hidden Markov model include the starting probability of the starting state. The transition probabilities from the starting state to insertion (I), deletion (D) and matching (M) are 1 / 3 respectively. The emission probabilities of the insertion parts are initialized to 1 / 4. The transition probabilities of all matching state positions are given the same initialization value, where TMM = 0.96, TMI = TMD = exp(-4), TDI= TID = 0, TII = TIM = TDM = TDD = 0.5, TM-end = 0.5.

[0064] 2. Spectral Hidden Markov Model for Single Sequences Each sequence can be calculated in parallel during processing without additional effort. However, for the directed acyclic pan-genome graph constructed from massive genomic data sets, since many nodes have multiple parent / child nodes, a dependency relationship determined by the topology order is formed. Therefore, the shingled spectral hidden Markov model is developed to execute parallel tasks that depend on the topology order. Figure 2 As shown in the figure, for a given training directed acyclic pan-genome graph, the traversal of the complementary order dependency is first carried out, and each path without branches in the graph is taken as a separate computing task, and the dependencies between different computing tasks are recorded. All computing tasks that meet the conditions are placed in a task queue. The idle thread extracts the task from the task queue, and after completion, it records and checks whether the subsequent tasks that depend on it can be executed. If all dependencies are satisfied, they are added to the end of the task queue.

[0065] 3. The classic profile hidden Markov model is developed for a single biological sequence and cannot be used directly on a directed acyclic graph. The main difference between a directed acyclic graph and a single sequence is that a node in the former can have multiple parent / child nodes, while each character in the latter has only one character before and after. The key to the present invention in meeting this challenge is to construct a single virtual parent / child node for a graph node with multiple parent / child nodes, thereby smoothly implementing the forward and backward calculations in the hidden Markov model training. Forward virtual node construction and probability summation calculation are as follows: Figure 3 After completing the sum of the virtual node probabilities, the forward probability calculation can be calculated using the conventional forward probability formula. The conventional forward probability calculation formula is as follows:

[0066]

[0067] Among them, f represents the forward probability, i represents the index of the node, i-1 represents the virtual parent node of the node, M, I, and D represent the matching state (Match), the insertion state (Insert), and the deletion state (Delete) respectively; k represents the index of the state, α represents the transition probability, Indicates that node i is in state The forward probability when X is M, I or D, Indicates that node i is in state The emission probability when Y is M or I, Indicates that the virtual parent node of node i is in state The forward probability when X is M, I or D, Indicates status Transfer to state The transition probability, X is M, I or D.

[0068] The backward virtual node is constructed as follows Figure 4 As shown, the probability sum formula is as follows:

[0069]

[0070] Where b represents the backward probability; i represents the index of the node; j represents the index of the child node; C represents the child node of the current node; n i represents the number of child nodes of node i; M, I, D represent the match state (Match), insert state (Insert) and delete state (Delete), respectively, k represents the state index; W represents the weight of the edge and the weight of the node, and e represents the emission probability. Indicates that the virtual child node of node i is in state The backward probability when X can be M, I or D, Indicates that the jth child of node i is in state The backward probability when X can be M, I or D, Indicates that the jth child of node i is in state The emission probability when X can be M or I, Represents the ratio of the weight of the edge from node i to its jth child node to the weight of node i.

[0071] After completing the construction of virtual nodes and probability summation, the backward probability calculation is calculated according to the conventional spectral hidden Markov model. The specific formula is as follows:

[0072]

[0073] Among them, b represents the backward probability, i represents the index of the node, i+1 represents the virtual child node of the node, M, I, and D represent the matching state (Match), the insert state (Insert), and the delete state (Delete), k represents the index of the state, and α represents the transition probability. For example Indicates that node i is in state The backward probability when , X can be M, I or D, Indicates that the virtual child node of node i is in state The forward probability when , X can be M, I or D. Indicates status Transfer to state The transition probability of , X can be M, I or D.

[0074] Compared with a single sequence set, a directed acyclic graph can achieve tens of thousands of effective compressions by merging identical fragments (the specific compression rate depends on the similarity / conservatism of the sequences, for example, the compression rate of 4 million high-quality COVID-19 genomes can reach 25,000 when the fragment length is 20). Correspondingly, for linear computational complexity comparison calculations, the amount of calculation is also reduced by approximately the multiple corresponding to the compression rate. For quadratic computational complexity, the upper limit of savings / acceleration can rise to the square of the compression rate. For cubic computational complexity, the upper limit of savings / acceleration continues to rise to the cube of the compression rate.

[0075] 4. In the classic spectral hidden Markov training technique, the computational time and space complexity of each sequence in each iterative calculation process is , where L is the sequence length, L PHMM is the number of matching states. This is an unbearable computing burden and memory requirement for sequences with longer lengths. The reason is that we have no idea about the relationships between individual sequences before completing the alignment, so we must consider the matching possibilities of each character in the sequence and each matching state in the spectral hidden Markov model. However, the directed acyclic pan-genome graph itself is a partially completed alignment. We first compare the alignment reference graph (see Figure 5 ) is defined as the reference path for comparison, and its length is defined as the number of matching states of the shingled hidden Markov model. In this way, the topological sequence coordinates in the reference path form a one-to-one correspondence with the matching states in the shingled hidden Markov model. The relative positions between different paths are constrained by the common upstream parent nodes or downstream child nodes of these paths. We use this existing information to implement a high-speed shingled training mode through two steps:

[0076] 1) Using the preprocessed comparison reference graph and reference path (see the comparison reference graph / path section below for details), determine the most likely matching range of non-reference path nodes on the reference path. The specific process is as follows: Figure 6 shown.

[0077] 2) Relative experience window between different reference path branches. Different path nodes that share topological sequence coordinates are likely to correspond to the same or similar spectral hidden Markov model matching states, so an experience window width that can maintain accuracy can be determined through fast calculation.

[0078] The sum of the maximum non-reference path node matching range and the experience window on the reference path is the shingle width W. That is, each node only considers matching with matching states within its shingle width range, and the probability of matching with matching states outside the shingle width range is set to zero. This reduces the computational space complexity from L PHMM L is reduced to L PHMMW, W depends on the similarity / conservation between sequences, and is much smaller than and has no direct relationship with the sequence length L. For the new coronavirus genome, we can obtain higher accuracy by taking W equal to 250, which saves about 120 times the computing time and memory compared to the sequence length of about 30,000.

[0079] 5. For conventional spectral hidden Markov models, Viterbi calculates the maximum probability state chain of reaching a specific character i when a single sequence is given a single sequence character (1,2,…,i-1). Since nodes in a directed acyclic pan-genome graph may have multiple parent nodes, similar to the forward and backward calculations in processing training, a virtual parent node (i-1) is constructed to implement the corresponding forward maximum probability calculation (such as Figure 6 As shown). In the Viterbi backtracking stage, when there are multiple parent nodes, only one parent node path will be determined as the maximum probability path in the initial backtracking, and the paths corresponding to other parent nodes need to be backtracked from the child nodes again. Another difficulty is that in the case of multiple child nodes, different genomes in the parent node can correspond to different hidden Markov model states. It is necessary to verify and exclude virtual paths and record the possible different states of different genomes in the same node.

[0080] Although the present invention has been described above with reference to the embodiments, various modifications may be made thereto and parts thereof may be replaced by equivalents without departing from the scope of the present invention. In particular, as long as there is no structural conflict, the various features in the embodiments disclosed in the present invention may be used in combination with each other in any manner, and the fact that these combinations are not exhaustively described in this specification is only for the sake of omitting space and saving resources. Therefore, the present invention is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A shingled hidden Markov sequence alignment method on a directed acyclic pan-genome graph, characterized in that: Here are the steps: S1. The directed acyclic pan-genome graph is subjected to pre-alignment pre-processing by reducing the fragment length, removing degenerate base nodes and low-weight nodes, and obtaining a Viterbi graph, a training graph and an alignment reference graph; S2, taking the length of the longest comparison reference path in the comparison reference graph as the number of matching states of the shingled spectral hidden Markov model, and training on the training graph after initializing the model parameters and the shingled width of each node to obtain the shingled spectral hidden Markov model, wherein the sum of the maximum non-reference path node matching range and the experience window on the reference path is the shingled width; S3. After obtaining the shingled spectral hidden Markov model, the Viterbi decoding algorithm based on virtual node probability calculation is used to calculate the most likely state path of each sequence in the shingled spectral hidden Markov model on the Viterbi graph, and the alignment relationship of each position in the sequence is determined according to the state path, thereby completing the multiple sequence alignment; In step S2, the shingled spectral hidden Markov model uses three different smoothing parameters: exp(-3), exp(-5) and exp(-7). The smoothing formula is as follows: ； Among them, for a selected set of n possible variables, p i is the probability of the ith possibility, p j is the probability of the jth possibility, c i is the smoothing parameter for the ith possibility, c j is the smoothing parameter for the jth possibility, is the updated probability.

2. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In step S1, the steps of reducing the length of the directed acyclic pan-genome graph segment are as follows: Decompose all starting nodes of the directed acyclic pan-genome graph into L1-L2+1 new nodes, and each node retains a sequence segment of length L2, where L1 represents a segment length of FTO-DAG; For all other nodes, only the last L2 characters are retained in the sequence fragments; Node merging based on topological sequence coordinates is performed in the graph, wherein no node merging is performed when the length of the last segment is reduced to 1 in the stepwise reduction of the segment length.

3. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In step S1, the step of removing degenerate base nodes from the directed acyclic pan-genome graph is as follows: all nodes and corresponding edges containing degenerate bases from the directed acyclic pan-genome graph are removed.

4. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In step S1, the steps of removing low-weight nodes in the directed acyclic pan-genome graph are as follows: First, all nodes below the global low weight threshold are removed from the directed acyclic pan-genome graph; Then, traverse from each tail node in reverse topology order, stop traversing at the first node whose weight is greater than the tail low weight threshold, and delete the tail low weight nodes that have been traversed; Among them, the global low weight threshold is set to 0.01 or 0.001, and the tail low weight threshold is set to 0.

01.

5. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In step S2, the parameters of the shingled spectral hidden Markov model include the starting probability of the starting state, the transition probabilities from the starting state to insert I, delete D and match M are 1 / 3 respectively, the emission probability of the inserted part is initialized to 1 / 4, and the transition probabilities of all matching state positions are given the same initialization value, where T MM = 0.96, T MI = T MD = exp(-4), T DI = T ID = 0, T II = T IM = T DM = T DD = 0.5, T M-end = 0.

5.

6. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In step S2, the shingled spectral hidden Markov model performs the following steps on the input directed acyclic pan-genome graph to be compared: First, we traverse the order-complement dependency, treat each path without branches in the directed acyclic pan-genome graph as a separate computing task, and record the dependencies between different computing tasks. All computing tasks that meet the conditions are placed in a task queue. The idle thread extracts the task from the task queue, and after completion, it records and checks whether the subsequent tasks that depend on it can be executed. If all dependencies are met, they are added to the tail of the task queue.

7. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In the shingled spectral hidden Markov model, for graph nodes with multiple parent / child nodes in the training data set, a single forward virtual parent / child node is constructed and probability sum calculation is performed, and a single backward virtual parent / child node is constructed and probability sum calculation is performed.

8. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In order to realize the parallel training of the shingled spectral hidden Markov model, each branchless path on the directed acyclic graph is taken as a computing task. The dependency between these tasks depends on their topological order. Tasks that are not dependent on each other are calculated simultaneously using different threads.

9. The method for shingled hidden Markov sequence alignment on a directed acyclic pan-genome graph according to claim 1, characterized in that: In the Viterbi decoding algorithm of the shingled spectral hidden Markov model, a virtual parent node i-1 is constructed for the graph nodes with multiple parent / child nodes in the data set to implement the corresponding forward maximum probability calculation, where i represents the specific character of a single sequence. In the backtracking stage of the Viterbi decoding algorithm, in the case of multiple parent nodes, only one parent node path is determined as the maximum probability path in the initial backtracking, and the paths corresponding to the other parent nodes are backtracked from the child nodes again. In the case of multiple child nodes, the virtual path is verified and excluded, and the possible different states of different genomes in the same node are recorded.

Images