A third generation sequence alignment method based on longest path search

Abstract

This invention provides a three-generation sequence alignment method based on longest path search. First, a hash index of the reference genome sequence is constructed. Then, each k-mer of the sequence to be aligned is extracted, and its position in the genome is found using the hash index. Each matching k-mer is treated as a node, and a directed acyclic graph of k-mer l-neighborhoods is constructed. Based on the position information of the matching k-mer in the sequence to be aligned, it can be determined whether nodes are connected by edges and their directions. Isolated nodes and small isolated networks are then filtered out. A dynamic scoring strategy is designed to determine the predecessor node of each node, select the highest score from the predecessor nodes, and record the score path. The longest path is selected. The sequence to be aligned and the reference genome can be divided into seed regions and non-seed regions. For non-seed regions, a traditional double sequence alignment method is used to obtain detailed base alignment results, which are finally merged with the seed regions to obtain the alignment results for the entire sequence.

Description

Technical Field

[0001] This invention relates to a method for processing long DNA sequences, which is mainly a three-generation sequence alignment method based on the longest path search. Background Technology

[0002] Third-generation sequences refer to sequences generated using third-generation DNA sequencing technology. Compared to second-generation sequencing sequences, third-generation sequences are longer, have higher error rates, and lower coverage. Third-generation sequence alignment involves comparing third-generation sequencing data with a reference genome to determine the similarity and location of the sequences. Sequence alignment is a crucial aspect of sequence analysis and forms the basis for subsequent bioinformatics mining. Most existing sequence alignment methods are designed for data generated by second-generation sequencing technology (shorter sequences, lower error rates) and are not suitable for processing third-generation sequence data. Therefore, new third-generation sequence alignment methods need to be developed. Currently, alignment methods for third-generation sequences mainly include hash-based seed search methods, text compression and indexing-based seed search methods, and seed search methods based on existing alignment tools. Each method has its advantages but also limitations. A common problem is that these methods have low alignment sensitivity and are highly sensitive to sequencing errors. Summary of the Invention

[0003] To overcome the shortcomings of existing methods, this invention provides a three-generation sequence alignment method based on the longest path search.

[0004] The purpose of this invention is to address the problems of long sequence lengths and high error rates in third-generation sequencing, as well as the low sensitivity and poor robustness of existing alignment methods to sequencing errors. This invention proposes a third-generation sequence alignment method based on the longest path search, which has high alignment sensitivity and coverage, and strong robustness to sequencing errors, thus providing effective technical support for the analysis of third-generation single-molecule sequencing data.

[0005] To achieve the above objectives, the basic idea of ​​this invention is as follows: First, a fast hash table is used to construct the k-mer position index information of the reference genome sequence. Then, each k-mer of the sequence to be aligned is extracted, and through hash function transformation, the positions of the k-mers of the sequence to be aligned in all parts of the genome can be easily and quickly found. Each matching k-mer is regarded as a network node, and a directed acyclic graph of k-mer l-neighborhood is constructed. Based on the position information of the matching k-mer in the sequence to be aligned, it can be determined whether two k-mers are connected by an edge and the direction of the edge. After the directed acyclic graph is established, isolated nodes and small isolated networks are filtered out. Then, a specific dynamic scoring strategy is designed to determine the predecessor node of each node, select the maximum score from the predecessor nodes, and record the scoring path. After obtaining the score value of each node, the node with the largest score value is selected, and then the longest path can be obtained by path backtracking. All nodes included in the longest path can divide the sequence to be aligned and the reference genome into seed regions and non-seed regions. For non-seed regions, traditional double sequence alignment methods are used for alignment to obtain detailed base alignment results. Finally, these results are merged with the seed regions to obtain the final alignment result of the entire sequence.

[0006] The three-generation sequence alignment method based on longest path search of this invention includes the following steps:

[0007] Step 1: Constructing a genome k-mer location library

[0008] A genome k-mer position library refers to a hash table that stores the positions of k-mer sub-fragments in the genome. First, all k-mer sub-fragments are extracted from the genome sequence. Then, a hash function is used to perform k-mer transformation, and the positions of the k-mers in the genome are stored in the hash table. The specific implementation process is as follows:

[0009] 1) Create a length of 4 based on the k-mer size (Γ). Γ A hash table, or array, is used to store the location of the corresponding k-mer in the genome;

[0010] 2) Extract all k-mers from the genome sequence. A k-mer is a sub-fragment of the genome sequence containing k bases. For a genome DNA sequence of length L, if the k-mer length is Γ, the total number of k-mers in the genome is L-Γ+1.

[0011] 3) Perform a hash transformation on the first k-mer of the genome sequence. Assume that this k-mer (of length Γ) can be represented as: w = c1, c2, ..., c Γ Its storage address in the hash table can be calculated using the following hash function:

[0012]

[0013]

[0014] The hash code of each k-mer can be seen as a transformation of a Γ-bit quaternary number. The hash value of each k-mer in the genome sequence is calculated by formula (1), which is the index position of the k-mer in the position library. Then, the position of this k-mer in the genome is stored in the array under the index.

[0015] 4) Repeat step 3) Calculate the hash value of all k-mers in the genome and store their positions in the genome. After storing all k-mer position information, the constructed genome k-mer position library is obtained.

[0016] Step 2: Matching and correcting k-mer representations

[0017] First, all k-mers of the sequence to be aligned are extracted; then, based on the genomic k-mer location library constructed in the previous step, the location of each k-mer in the genome of the sequence to be aligned is found; then, the location information of each match is corrected so that the distribution trend of the reverse matching nodes is consistent with the trend of the forward matching nodes; the specific implementation steps are as follows:

[0018] 1) For the sequence r to be aligned, extract all k-mers of r, and then calculate the hash value using formula (1) to find the location information of the k-mer in the genome. Each k-mer of sequence r and its location in the genome can be represented by a quadruple:

[0019] match i (pr,pg,dr,dg)(3)

[0020] In the formula, match i (pr) represents the position of the i-th match in sequence r, match i (pg) represents the position of the i-th match in the reference genome. i (dr) and match i (dg) represents the direction of the i-th match in sequence r and in the reference genome (1 indicates forward direction, 0 indicates reverse direction);

[0021] 2) After finding the quadruple information of each k-mer in sequence r using formula (3), the position information of the k-mer in reverse matching is corrected using the following formula, and converted into triple information. This makes the changing trend of the k-mer in reverse matching consistent with the changing trend of the k-mer in forward matching, and linearizes the k-mer in correct matching. The correction formula is:

[0022]

[0023] It can be observed from the above formula that m i and match i The main difference is the representation of the position and direction in the sequence to be compared; if match i (dr) and match i (dg) are in the same direction, then m i (r) = '+', m i (pr) = match i (pr), otherwise m i (r) = '-', m i (pr) = len(r) – match i (pr), where len(r) represents the length of the sequence r; m i (r) represents the consistency of the matching direction, '+' indicates that the current k-mer matches the positive strand of the reference genome, '-' indicates that it matches the reverse strand; finally, the corrected set of k-mers is represented as: M = {m1,..., m i ,..., m Nm}, where Nm represents the total number of k-mers;

[0024] Step 3: Construct a k-mer directed acyclic graph

[0025] Assume that all k-mers in M can be used to construct a directed unweighted network, where each k-mer can be used as a node, and an edge is determined to be connected between each pair of k-mers according to their distance; for each node v i (v i = m i ), calculate its l-neighborhood nodes, and add an edge from v i to its neighborhood nodes; for node v i , its l-neighborhood nodes are defined as the set of nodes that satisfy the following formula: v j (pr) - v i (pr) < l. The parameter l represents the maximum allowed distance between two nodes. Generally, l is the length of the sequence to be compared.

[0026] Step 4: Selection of the longest path

[0027] For the directed graph constructed above, first filter out single nodes and smaller isolated subgraphs (subgraphs with fewer than 8 nodes), which can avoid the path search space. For the remaining nodes, select the longest path containing the largest number of nodes as the best alignment skeleton; the selection of the longest path first scores each node, and the formula is:

[0028]

[0029] In the formula, score(v i ) represents the score of each node, pre(v i ) indicates pointing to v i The set of nodes (or predecessor nodes); representing the selection of the maximum score from the predecessor nodes of each node and recording the score path; considering the possible sequencing errors or structural variations in the sequence, the top N (the default value of N is 10) longest non-overlapping paths are selected as the final alignment backbone;

[0030] Step 5: Final Comparison Results

[0031] Each path contains nodes that can segment the sequence to be aligned and the reference genome sequence into paired sub-fragments. Then, the traditional double sequence alignment method is used to obtain detailed alignment results of the sub-fragments, which are finally merged with k-mer to obtain the alignment results of the entire sequence.

[0032] Preferably, step 2, by correcting the information of each matching k-mer, can easily make the position change trend of the reverse matching k-mer consistent with the position change trend of the forward matching k-mer. This is beneficial for constructing a directed acyclic graph of k-mers and scoring them according to a unified formula, and facilitates the search for the longest forward k-mer path and the reverse matching path.

[0033] Preferably, step 3 determines the connection edges and directions of each node through l-neighborhood. This can maximize the continuity of the path, search for paths in regions with high sequencing errors or structural variations, avoid the generation of fully connected networks, reduce the time complexity of the longest path search, and improve the overall efficiency of the method.

[0034] The present invention has the following beneficial effects:

[0035] By treating the matched k-mers as nodes and correcting the positions of the reverse-matched k-mers, and then determining the alignment frame by constructing a directed acyclic graph and searching for the longest path, the method improves the sensitivity of sequence localization and enhances its robustness to sequencing errors. Attached Figure Description

[0036] Figure 1 This is a flowchart of the three-generation sequence alignment method based on the longest path search of the present invention.

[0037] Figure 2 This is a schematic diagram of the reverse matching k-mer correction;

[0038] Figure 3 This is a schematic diagram showing the alignment results of different methods under different sequencing errors;

[0039] Figure 4This is a diagram illustrating the comparison results of different methods. Detailed Implementation

[0040] On a server with an Intel Xeon Gold 5218R CPU@2.10GHz and 256GB of RAM, based on a CentOS 7.5 Linux platform, eight sets of simulated data with different error rates and three sets of real sequencing datasets from different sequencing platforms were selected for comparative simulation experiments. The simulated datasets were generated by simulation software, and the real datasets were generated by PacBio, MinONT, and Illumina sequencing platforms.

[0041] The present invention will be further described below with reference to the accompanying drawings and embodiments. Figure 1 This is a flowchart of the three-generation sequence alignment method based on the longest path search of the present invention; Figure 2 This is a schematic diagram of reverse matching k-mer correction. The left diagram shows the original matching position information, and the right diagram shows the corrected matching position information. Figure 3 This is a schematic diagram of alignment results for different methods under different sequencing errors; the horizontal axis represents sequencing error, the vertical axis of the upper graph represents the percentage of correctly aligned sequences, and the vertical axis of the lower graph represents the proportion of correctly aligned bases. Figure 4 This is a schematic diagram showing the alignment results of different methods. The horizontal axis represents the number of aligned sequences, and the vertical axis represents different reference genomes.

[0042] If the k-mer length is set to 11, the specific alignment steps are illustrated using the first sequence r (2824bp in length) from the E. coli sequence dataset.

[0043] Step 1: Constructing a genome k-mer location library

[0044] 1) Since the k-mer length is set to 11, the creation length is 4. 11 An array of (4194304);

[0045] 2) Extract all k-mers from genome g. Since g is 4,681,865 bp long and k-mers are 11, the total number of k-mers is 4,681,855.

[0046] 3) Calculate the hash value of the first k-mer in genome g. The first k-mer sequence is CAAGCCAGCCA, and the hash value can be obtained as 1086612 using formula (1). Since this k-mer is the first k-mer in the genome, store 0 in the hash table at the 1086612th position.

[0047] 4) Repeat step 3) until all k-mer hash values ​​of the genome are calculated and their positions in the genome are stored in a hash table, which is the constructed genome k-mer position library.

[0048] Step 2: Find matching k-mer and correct positional information

[0049] Extract the first k-mer of sequence r, w1 = ATCACCACAAG, its reverse sequence is CTTGTGGTGAT. According to formula (1), its forward and reverse hash values ​​are obtained. By searching the corresponding position elements of the genome k-mer position library, it is found that the reverse k-mer (hash value 2079459) appears once in genome g, at position 445058. According to formula (3), the quadruple information of w1 can be constructed: match i (1,445058,-,+), according to formula (4), we correct it to obtain m. i =(2823,445058,-).

[0050] Step 3: Calculate each m i The neighboring nodes are identified and the edges are determined. Single nodes and small isolated subgraphs (subgraphs with fewer than 8 nodes) are filtered out. For the remaining nodes, the longest path containing the most nodes is selected as the best alignment skeleton according to formula (5).

[0051] Step 4: Use traditional double sequence alignment methods, such as the Needleman-Wunsch method, to align the unaligned regions, and then merge them with the seed region to obtain the final alignment result. Figure 1 As shown.

[0052] Table 1 shows the results of different methods on three sets of real data. Four indicators were used for comparison: the number of aligned sequences, the number of aligned bases, the total alignment score, and the number of matching bases. The number of aligned sequences and the number of aligned bases reflect alignment sensitivity, while the total alignment score and the number of matching bases reflect alignment quality. As shown in Table 1, the longest path search-based third-generation sequence alignment method achieved high values ​​for each indicator, indicating that it has high alignment sensitivity and alignment quality.

[0053] Table 1. Sequence alignment results of three generations based on longest path search, compared with other methods.

[0054]

[0055]

[0056] The above results demonstrate that the third-generation sequence alignment method based on longest path search can perform sequence alignment on third-generation sequencing data with long sequences and high sequencing error rates, yielding more aligned sequences and bases while maintaining good alignment quality. It is suitable for aligning various types of genome third-generation sequencing data and has significant potential application value.

Claims

1. A method for aligning three sequences based on a longest path search, characterized in that, Includes the following steps: Step 1: Constructing a genome k-mer location library A genome k-mer position library refers to a hash table that stores the positions of k-mer sub-fragments in the genome. First, all k-mer sub-fragments are extracted from the genome sequence. Then, a hash function is used to perform k-mer transformation, and the positions of the k-mers in the genome are stored in the hash table. The specific implementation process is as follows: 1) Create a hash table, i.e. an array, of length 4 Г for each k-mer size Г to store the positions of the corresponding k-mers in the genome; 2) Extract all k-mers from the genome sequence. A k-mer is a sub-fragment of the genome sequence containing k bases. For a genome DNA sequence of length L, if the k-mer length is Г, the total number of k-mers in the genome is L-Г+1. 3) Hash conversion is performed on the first k-mer of the genome sequence, assuming that the k-mer with a length of Г is represented as: The storage address of the k-mer in the hash table is calculated by the following hash function: （1） In formula 4 Г-γ It is the first in k-mer Bases at each position ( The weight of ) It is an indexing function, defined as: （2） The hash code of each k-mer is regarded as A conversion of the binary number to the quaternary number, the hash value of each k-mer of the genome sequence is calculated by formula (1), which is the index position of the k-mer in the position library, and then the position of the k-mer in the genome is stored in the array under the index. 4) Repeat step 3) Calculate the hash value of all k-mers in the genome and store their positions in the genome. After storing all k-mer position information, the constructed genome k-mer position library is obtained. Step 2: Matching and correcting k-mer representations First, all k-mers of the sequence to be aligned are extracted; then, based on the genomic k-mer location library constructed in the previous step, the location of each k-mer in the genome of the sequence to be aligned is found; then, the location information of each match is corrected so that the distribution trend of the reverse matching nodes is consistent with the distribution trend of the forward matching nodes; the specific implementation steps are as follows: 1) For the sequence r to be aligned, extract all k-mers of r, and then calculate the hash value using formula (1), that is, find the location information of the k-mer in the genome. Each k-mer of sequence r and its location in the genome is represented by a quadruple: match i (pr, pg, dr, dg) (3) In the formula, match i (pr) represents the position of the i-th match in sequence r, match i (pg) represents the position of the i-th match in the reference genome, match i (dr) and match i (dg) represents the orientation of the i-th match in sequence r and in the reference genome, respectively, with 1 indicating forward and 0 indicating reverse; 2) After finding the quadruple information of each k-mer in sequence r using formula (3), the position information of the k-mer in reverse matching is corrected using the following formula, which converts it into triple information, and makes the changing trend of the k-mer in reverse matching consistent with the changing trend of the k-mer in forward matching; so as to linearize the k-mer in correct matching; the correction formula is: （4） Based on the above formula, m was observed... i With match i The main difference lies in the position and orientation representation in the sequences to be compared; if match i (dr) and match i If the directions (dg) are the same, then m i (r) = '+', m i (pr) = match i (pr), otherwise m i (r) ='-',m i (pr) = len(r) – match i (pr), where len(r) represents the length of sequence r; m i (r) indicates consistency in matching direction, '+' indicates that the current k-mer matches the positive strand of the reference genome, and '-' indicates that it matches the negative strand; finally, the corrected k-mer set is represented as: , where Nm represents the total number of k-mer matches; Step 3: Construct a k-mer directed acyclic graph Suppose that all k-mers in M ​​are used to construct a directed unweighted network, where each k-mer is a node, and each pair of k-mers is connected by an edge based on their distance; for each node v i v i = m i Calculate its l-neighborhood nodes and add v i Edges pointing to its neighboring nodes; for node v i Its l-neighborhood nodes are defined as the set of nodes that satisfy the following formula: v j (pr) - v i (pr) < l; the parameter l represents the maximum allowed distance between two nodes, and in general l is the length of the sequence to be compared; Step 4: Selecting the longest path For the directed graph constructed above, isolated subgraphs with fewer than 8 nodes are first filtered out to avoid path search space. For the remaining nodes, the longest path containing the most nodes is selected as the best alignment skeleton. The longest path selection is done by scoring each node using the following formula: （5） In the formula, score(v i ) represents the score of each node, pre(v i ) indicates pointing to v i The set of nodes is called the predecessor node; it means selecting the maximum score from the predecessor nodes of each node and recording the score path; considering the sequencing errors or structural variations in the sequence, the first N longest non-overlapping paths are selected as the final alignment backbone, and the default value of N is 10. Step 5: Final Comparison Results Each path contains nodes that segment the sequence to be aligned and the reference genome sequence into paired sub-fragments. Then, the traditional double sequence alignment method is used to obtain detailed alignment results of the sub-fragments, which are finally merged with k-mer to obtain the alignment results of the entire sequence.

2. The method of claim 1, wherein the method is based on a longest path search. Step 2 corrects the information of each matching k-mer, making it easier to ensure that the position change trend of the reverse matching k-mer is consistent with the position change trend of the forward matching k-mer. This facilitates the construction and scoring of the k-mer directed acyclic graph according to a unified formula, and makes it easier to find the longest k-mer forward path and the reverse matching path.

3. The method of claim 1, wherein the method is based on a longest path search. The step 3 determines the edges and directions of each node through the l-neighborhood, guarantees the continuity of the path, searches the path of the region with high sequencing error or the structural variation region, avoids the generation of the full connection network, reduces the time complexity of the longest path search, and improves the operation efficiency of the whole method.

Images