CN103699818B - Two-way side extended method based on the elongated kmer inquiries of the two-way De Bruijns of multistep - Google Patents

Abstract

The present invention relates to the technical field of gene sequencing, and provides a vertex expansion method for variable-length kmer query based on a multi-step bidirectional De Bruijn graph, comprising: step A: reading a sequencing data source file, and constructing a multi-step bidirectional De Bruijn graph; B: performing statistics on intersecting bidirectional edges in the multi-step bidirectional De Bruijn graph; Step C: expanding bidirectional edges based on variable-length kmer queries in the multi-step bidirectional De Bruijn graph. The present invention only constructs variable-length kmer combinations for existing bifurcated edges on the De Bruijn graph, and queries its occurrence times in the input sequence, and then selects the merging of bidirectional edges according to the times; selects the optimal merging of bifurcated edges , to minimize the possibility of erroneous merging of bidirectional edges; it can significantly increase the length of contigs, and can also minimize the quality loss of contigs.

Description

Bidirectional edge expansion for variable-length kmer queries based on multi-step bidirectional De Bruijn graphs method

【Technical field】

The invention relates to the technical field of gene sequencing, in particular to a bidirectional edge extension method of variable-length kmer query based on a multi-step bidirectional De Bruijn graph.

【Background technique】

Gene sequence analysis takes algorithms and mathematical models as the core, and the research content involves many aspects, mainly including: storage and acquisition of genetic data, sequence comparison, sequencing and splicing, gene prediction, biological evolution and phylogenetic analysis, protein structure prediction, RNA structure prediction, molecular design and drug design, metabolic network analysis, gene chip, DNA computing, etc. Now the close combination of biotechnology and computer information processing technology has accelerated the speed of processing biological information data, making the interpretation of biological meaning as accurate as possible in the shortest possible time, and accelerating the development of bioinformatics. At present, biological information processing has become one of the great challenges faced by the current information technology field.

Gene sequence analysis is the analysis of massive gene sequence data to extract and mine new biological information knowledge. Among them, it involves machine learning, pattern recognition, book analysis and mining, combinatorial mathematics, random models, strings, graph algorithms, distributed computing, high-performance computing, parallel computing and other knowledge in computer technology. Among them, the study of whole genome is one of the cores of current bioinformatics research.

Gene is the most basic genetic code of human beings, representing the life information of each person. There are subtle differences in genetic loci in the gene sequence, and these genetic code polymorphisms are closely related to human health, pathogenic mechanism, and medical treatment. Among them, DNA sequencing is one of the basic contents that need to be completed to study the whole genome sequence.

Since the advent of Sanger sequencing technology in 1977, after more than 30 years of development, DNA sequencing technology has developed by leaps and bounds. The second-generation sequencing technology characterized by high-throughput and short sequences has gradually occupied the market. Three generations of sequencing technologies have gradually emerged, and they each have different advantages in sequencing characteristics. The data extraction and analysis software of the traditional gene sequencing method has been relatively perfect after nearly 10 years of research and development. However, the development of sequencing technology has brought about changes in sequencing data, making the existing data processing software unable to meet the needs of current biomedical research.

The application of next-generation high-throughput sequencing methods can complete the determination of the entire genome data in a short time. The rapid development of high-throughput sequencing methods also poses challenges to the analysis and processing methods of the acquired genetic data. In this currently hot research field, there is an urgent need to develop a bioinformatics platform that can meet the massive data processing requirements of high-throughput sequencing technologies. In the face of the Personal Genome Project and the prospect of personalized medicine in the future, high-efficiency and low-cost sequencing technology has become an inevitable trend. At the same time, complete sequencing solutions such as simplified and efficient one-stop complete bioinformatics data analysis platforms are also extremely important and indispensable development directions.

However, although the next-generation high-throughput sequencing method has a high sequencing throughput, it will introduce sequencing errors. At the same time, the sequencing errors of the sequencing sample itself, uneven sequencing, SNP, and the repeated sequence Repeat of the genome itself, and these sequencing errors, SNP , Repeated sequences will introduce some wrong bidirectional edges or vertices in the multi-step bidirectional De Bruijn graph constructed during genome assembly, making many bidirectional edges unable to continue to expand. At the same time, due to the fixed kmer length, the effective information of the sequencing sequence is lost, and it is impossible to decouple all repetitive sequences whose length exceeds the kmer length. The above conditions make it impossible for the De Bruijn graph to continue to shrink, and the contig cannot be expanded, which ultimately makes the length and quality of the contig very low.

The assembly of short gene fragments generated by next-generation high-throughput sequencing methods leads to a large number of sequencing errors, which greatly increases the computational load of assembly algorithms. A large number of sequencing errors increase the assembly error rate and seriously affect the assembly results. Whether this problem can be effectively solved becomes the key to evaluating the pros and cons of an assembly algorithm.

At present, the strategy of assembly algorithm is mainly divided into two categories, one is the algorithm based on Overlap-Layout-Consensus (OLC), and the other is the algorithm based on De Bruijn graph. Among them, the software developed based on the OLC assembly algorithm, such as SSAKE, VCAKE, SHARCGS, etc., has more advantages in the assembly of long gene sequences, but it is not fully suitable for the new generation of short sequence assembly. Different from the OLC assembly algorithm, the De Bruijn algorithm no longer organizes data in units of reads, but assembles data in units of k-mers. Its advantages mainly include the following aspects: First, sequence assembly is performed in units of k-mers , does not affect the quality of nodes, and reduces the amount of redundant data. Secondly, the repeated area only appears once in the figure, which is easy to identify, can avoid wrong assembly, and reduce the error rate. Finally, a strategy of mapping overlapping regions to the same arc is adopted, which simplifies the search path. At present, many short sequence assembly algorithms use this framework, such as Velvet, IDBA, SOAPdenovo, ABySS, etc.

Among them, Velvet effectively utilizes the De Bruijn graph to realize efficient short sequence assembly. Velvet constructs a De Bruijn graph with k-mer as the basic unit, uses the structure of the graph, combines the corresponding sequence features, simplifies the structure of the graph, and finally finds an optimal path to complete the assembly process. Velvet focuses on three structures generated by erroneous data, namely tip, bubble, and erroneous connection. According to the length principle and the minority principle, it removes all the lengths less than 2k; uses the depth-first search strategy in the Tour Bus algorithm to merge bubbles, and finally uses the coverage threshold method to remove erroneous connections. This method also makes full use of the paired-end double-ended information to further solve the repeat problem and optimize the assembly effect. Velvet makes full use of the structural properties of graphs, simplifies data redundancy, and has a great improvement in speed compared with previous algorithms. Although it does not correct the sequence in the preprocessing stage, its error prevention mechanism largely makes up for this defect. This makes it better for use in the assembly of large genome sequences.

IDBA is also based on the De Bruijn graph, which realizes simple and efficient short sequence assembly. IDBA uses k-mer as the basic unit. Unlike the past, it uses a variable k value range (Kmin-Kmax) instead of using a fixed k value to obtain the length of k-mers. Since gene assembly is based on k-mers, many overlapping units are usually formed, which makes the assembly face the problems of wrong position assembly, missing vertices and low coverage. The correct selection of the value of k becomes a key factor in assembly. The generation of some wrong reads also led to a large number of branching. The smaller the K value, the more serious the branching problem, and the larger the k value, the fewer repeat areas appear, which directly affects the quality of assembly. IDBA uses an unfixed k value for assembly, which can solve the branching problem well, thereby improving the quality of assembly. In addition, IDBA reduces the memory usage of IDBA significantly by deleting the wrong k-mers with low coverage, and also improves the processing speed of IDBA.

SOAPdenovo can efficiently assemble hundreds of millions of reads with high quality. SOAPdenovo inherits the advantages of OLC algorithm and De Bruijn graph algorithm, which greatly improves its assembly quality. SOAP reduces the occurrence of error sequences by means of filtering and error correction by presetting the k-mer threshold. At the same time, it successfully handles bubbles by referring to the method of Velvet software, which increases its average coverage. In addition, SOAPdenovo uses double-ended information for overlapping region matching, and merges reads to generate contig fragments to generate a contig-based graph structure. Thus, SOAPdenovo greatly simplifies the complexity of the contig graph.

ABySS introduces the idea of parallel computing, builds a linux cluster, and builds a distributed De Bruijn graph structure on the cluster to store data distributedly on each node. It uses the MPI communication mechanism to complete the mutual communication between nodes. From building graphs, error correction processing to subsequent fixed-point fusion, and finally completing the reproduction of the entire genome sequence, it has a great advantage in terms of running time and memory consumption, and its error rate is extremely low. In terms of performance, especially in a single machine in a cluster The use of memory has been greatly improved, and it is being used more and more widely.

The existing main sequence assembly software, such as SOAPdenovo, Velvet, ABySS, Ray, etc., are all based on the kmer of a given length for De Bruijn construction and then contraction. The only way to optimize it is to choose the best kmer length. This fixed-length kmer-based assembly strategy cannot be decoupled for all repeat sequences with a length of approximately kmer length. Although IDBA can iteratively shrink the De Bruijn graph for various kmer lengths, it needs to decompose, store and calculate all the sequences for each kmer length, and this strategy will consume huge memory and computing time.

In view of this, it is an urgent problem to be solved in this technical field to overcome the defects in the prior art.

【Content of invention】

The technical problem to be solved by the present invention is to provide a bidirectional edge extension method of variable length kmer query based on multi-step bidirectional De Bruijn graph.

The present invention adopts following technical scheme:

A bidirectional edge extension method for variable-length kmer queries based on multi-step bidirectional De Bruijn graphs, including:

Step A: Read the sequencing data source file and construct a multi-step bidirectional De Bruijn graph;

Step B: performing statistics on intersecting bidirectional edges in the multi-step bidirectional De Bruijn graph;

Step C: Bidirectional edge expansion based on variable length kmer queries in the multi-step bidirectional De Bruijn graph.

Further, in the step B, assemble into a longer kmer, and query the occurrence times of the correspondingly longer kmer, and use the size of the occurrence times as the weight to select and merge the direction of the bifurcation edge.

Further, in the step C, after converting all possible bifurcations of each forked bidirectional edge into a longer kmer, according to the number of occurrences, the fork path with the highest number of occurrences is used as the final The bidirectional edges on the corresponding forked paths are deleted together after merging.

Further, the step B further includes:

Step B1: traverse each bidirectional edge e in the De Bruijn graph;

Step B2: count the maximum length of the sequencing sequences of all input files, and assign it to Lmax;

Step B3: If the length of the two-way side is greater than the length of Lmax-kmer, return to step B1, otherwise execute step B4;

Step B4: The starting vertex of the bidirectional edge e is u, the starting direction is Direc_u, the ending vertex is v, and the ending direction is Direc_v;

Step B5: Take the last k+1 characters of all bidirectional edges that enter vertex u and end in Direc_u and put them into the string array m, and take the last character of all bidirectional edges that leave vertex v and start in Direc_v put into string array n;

Step B6: Combining the left and right sides of the bidirectional sides with the strings in the m and n arrays to construct a variable-length kmer array, wherein each element in the variable-length kmer array is variable-length kmer=m+e+n.

Further, the step C further includes:

Step C1: Open the sequence file and read each sequence one by one;

Step C2: read the sequence (including the match of its reverse sequence) of each matching sequence of the variable-length kmer set, and count each variable-length kmer;

Step C3: traverse each bidirectional edge e in the De Bruijn graph;

Step C4: Count the maximum length of the sequencing sequences of all input files, and assign it to Lmax;

Step C5: If the length of the two-way side is greater than the length of Lmax-kmer, return to step C3, otherwise execute step C6;

Step C6: The starting vertex of the bidirectional edge e is u, the starting direction is Direc_u, the ending vertex is v, and the ending direction is Direc_v;

Step C7: Take all the edges that enter the vertex u and end in the direction of Direc_u, take the last k+1 characters of the last two-way edge e_i and put them into the string array m, and take all the two-way edges e_j that leave the vertex v and start in the direction of Direc_v The last character is put into the string array n;

Step C8: Combine the left and right sides of the bidirectional side e with the strings in the m and n arrays to construct a variable-length kmer array, wherein each element in the variable-length kmer array is variable-length kmer=m+e+n. Query these variable-length kmers in the counted variable-length kmers, and take out the one with the largest count value, and merge the corresponding three sides e_i, e, and e_j, and then delete e_i, e, and e_j.

Further, in the step C2, the sequence read into each match of the variable-length kmer set includes a match to the reverse sequence.

Further, said step A further includes:

Compress storage steps, specifically

A11. Read a sequence s;

A12, using a sliding window to cut the sequence s into multiple fragments t;

A13, for each fragment t, use the nucleic acid coding table to encode, and express as a 64-bit integer a;

A14. Reverse the fragment t, use the symmetrical complementation table to complement the reversed fragment to obtain the complementary fragment v, and use the nucleic acid coding table in step A13 to encode the complementary fragment again, and express it as a 64-bit integer b;

A15, take the maximum number of integer a and integer b, as the mark number of k molecules of fragment t and complementary fragment v;

A16. Steps A11-A15 are repeated until all sequences are completed;

and De Bruijn graph construction steps, specifically as

A21. Read a sequence s;

A22, the sequence s is cut into a plurality of fragments t by a sliding window, a fragment t is selected as cur, and the marks of the preceding and following fragments are respectively pre and lat;

A23, if the encoding of t is less than its complementary fragment encoding, then exchange the values of pre and lat;

A24, the corresponding bit position in the forward position mapping table of cur is 1 to indicate the edge pointing to pre;

A25. The corresponding bit position in the reverse position mapping table of cur is 1 to indicate the edge pointing to lat;

A26. Repeat steps A22-A25 to process other fragments t of the sequence s until all fragments t of the sequence s are completed, and then execute step A27;

A27, read a new sequence s, repeat steps A22-A26; until all sequences are processed, execute step A28;

A28. Complete the construction of the bidirectional multi-step De Bruijn graph.

Further, the sliding window in the steps A12 and A22 is a sliding window with a length of k, where 0<k<32 and k is an odd number;

The nucleic acid coding table in the step A13 is {A:00, C:01, G:10, T:11};

The symmetrical complementary table in the step A14 is {A->T, C->G, G->C, T->A};

The step A14 specifically is to reverse the character string of the segment t, use a symmetric complementary table to change each character in the reversed character string into its complementary character, obtain the character string v of the complementary character, and use step A13 again The nucleic acid encoding table in encodes the string v and expresses it as a 64-bit integer b.

Further, in the step A22, if there is no segment before or after the segment t, the value of pre or lat is assigned empty or NULL;

In step A24, the forward position mapping table is {A:0, C:1, G:2, T:3}, and the position query character is the last character of pre;

In step A25, the reverse position mapping table is {A:4, C:5, G:6, T:7}, and the position query character is the complementary character of the first character of lat.

Compared with prior art, the beneficial effect of the present invention is:

(1) Construct a variable-length kmer combination only for the existing bifurcation edges on the De Bruijn graph, and query the number of occurrences in the input sequence, and then select the merge of two-way edges according to the number; while IDBA iterates all kmer length, it is necessary to construct all possible kmers of each kmer length, and then shrink the De Bruijn graph, which method will lead to greater memory consumption and calculation time consumption;

(2) Select the optimal merging of bifurcated edges to minimize the possibility of erroneous merging of bidirectional edges;

(3) The length of contigs can be significantly increased, and the loss of contig quality can also be minimized; compared with other existing methods that increase the contig length and sacrifice contig quality, the present invention has a certain degree of control and improvement.

【Description of drawings】

Fig. 1 is the flow chart of the bidirectional edge extension method of the variable length kmer query based on the multi-step bidirectional De Bruijn graph according to the embodiment of the present invention;

Fig. 2 is the flow chart of the compression storage step in Fig. 1 step A;

Fig. 3 is a flowchart of the construction steps of De Bruijn diagram in Fig. 1 step A;

Fig. 4 is the flow chart that step B carries out statistics on crossing bidirectional edges in the multi-step bidirectional De Bruijn graph;

FIG. 5 is a flow chart of step C of bidirectional edge expansion based on variable-length kmer queries in the multi-step bidirectional De Bruijn graph.

【detailed description】

In order to make the object, technical solution and advantages of the present invention more clear, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention.

In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not constitute a conflict with each other.

The purpose of the present invention is to design an expansion method based on variable-length kmer query crossing two-way edges, which will make the De Bruijn graph continue to shrink and contigs continue to expand, and will not introduce errors at the same time, resulting in the decline of contig quality and accuracy.

The present invention provides a kind of bidirectional edge expansion method based on the variable length kmer query of multi-step bidirectional De Bruijn graph, as shown in Figure 1, the method comprises:

Step A: Read the sequencing data source file and construct a multi-step bidirectional De Bruijn graph;

Step B: counting the intersecting bidirectional edges in the multi-step bidirectional De Bruijn graph;

Step C: Bidirectional edge expansion based on variable-length kmer queries in multi-step bidirectional De Bruijn graphs.

Among them, step A can be specifically realized in the following ways:

In the compressed storage step, the required raw data includes the FASTA format files generated by the first-generation, second-generation and next-generation sequencing instruments. The sequences in the FASTA files are cut into k molecules one by one and compressed and stored in binary code as a 64 Long integer k number of bits in the numerator.

As shown in Figure 2, specifically

A11, read a sequence s; wherein, the sequence s is taken from the FASTA format file;

A12, using a sliding window to cut the sequence s into multiple fragments t;

A13, for each fragment t, use the nucleic acid coding table to encode, and express as a 64-bit integer a;

A14, reverse the fragment t, use the symmetrical complementation table to complement the reversed fragment to obtain a complementary fragment, and use the nucleic acid coding table in step A13 to encode the complementary fragment again, and express it as a 64-bit integer b ;

A15, take the maximum number of integer a and integer b, as the mark number of k molecules of fragment t and complementary fragment v;

A16. Repeat steps A11-A15 until all sequences are completed.

Through the above steps, the kmer in the two traditional De Brujin diagrams is converted into a 64-bit k-molecule sign number for storage. This step can store the two compressed kmers in other software such as velvet, IDBA, and SOAPdenovo as a compressed number of k molecules, and after obtaining the number of k molecules, the length of the k molecule can also be calculated as k The fragment t and its complementary fragment v.

The sliding window in steps A12, A22 is a sliding window whose length is k, wherein 0<k<32 and k is an odd number; the nucleic acid coding table in step A13 is {A:00, C:01, G:10, T: 11}; the symmetric complementary table in step A14 is {A->T, C->G, G->C, T->A}; step A14 is specifically to reverse the character string of segment t, using symmetric The complementary table changes each character in the reversed character string into its complementary character to obtain the character string v of the complementary character, and uses the nucleic acid encoding table in step A13 to encode the character string v again, and expresses it as a 64-bit integer b.

and De Bruijn graph construction steps, 1, use the mark number of calculating k molecule in the above-mentioned compressed storage step, 2, use each fragment and the extended characters of its adjacent fragments as this k molecule and its adjacent fragments 3. Store the initialized k-molecule data structure into the hash_map with the key value of the k-molecule number as the key value.

As shown in Figure 3, specifically

A21. Read a sequence s;

A22, the sequence s is cut into a plurality of fragments t by a sliding window, a fragment t is selected as cur, and the marks of the preceding and following fragments are respectively pre and lat;

A23, if the coding of t is smaller than its complementary fragment coding, then exchange the values of pre and lat;

A24. The corresponding bit position in the forward position mapping table of cur is 1 to indicate the edge pointing to pre;

A25. The corresponding bit position in the reverse position mapping table of cur is 1 to indicate the edge pointing to lat;

A26. Steps A22-A25 are repeated to process other segments t of the sequence s until all segments t of the sequence s are completed, and step S27 is executed;

A27, read a new sequence s, repeat steps A22-A26; until all sequences are processed, execute step A28;

A28. Complete the construction of the bidirectional multi-step De Bruijn graph.

In step A22, if the fragment t has no previous or subsequent fragments, the value of pre or lat is assigned empty or NULL; the forward position mapping table in step A24 is {A:0, C:1, G:2, T:3 }, the position query character is the last character of pre; the reverse position mapping table is {A:4, C:5, G:6, T:7} in step A25, and the position query character is the first character of lat Complementary characters.

In step B, assemble into a longer kmer for all possible cases where the bidirectional edges of all left and right vertices in the multi-step bidirectional De Bruijn graph are merged with the bifurcated edges of the bidirectional edges, and query the corresponding The number of occurrences of the longer kmer is used as the weight to select and merge the direction of the bifurcation edge.

In this method, the data structure of each vertex in the multi-step bidirectional De Bruijn graph is set as:

As shown in Figure 4, step B further includes:

Step B1: traverse each bidirectional edge e in the De Bruijn graph;

Step B2: count the maximum length of the sequencing sequences of all input files, and assign it to Lmax;

Step B3: If the length of the two-way side is greater than the length of Lmax-kmer, return to step B1, otherwise execute step B4;

Step B4: The starting vertex of the bidirectional edge e is u, the starting direction is Direc_u, the ending vertex is v, and the ending direction is Direc_v;

Step B5: Take the last k+1 characters of all bidirectional edges that enter vertex u and end in Direc_u and put them into the string array m, and take the last character of all bidirectional edges that leave vertex v and start in Direc_v put into string array n;

Step B6: Combining the left and right sides of the bidirectional sides with the strings in the m and n arrays to construct a variable-length kmer array, wherein each element in the variable-length kmer array is variable-length kmer=m+e+n. When there are still two-way edges to be processed, return to step B1; when there are no two-way edges to be processed, end.

In step C, after converting all possible bifurcations of each bidirectional edge into a longer kmer, the fork path with the highest number of occurrences is used as the final fork merging scheme according to the number of occurrences , and the bidirectional edges on the corresponding bifurcated paths are deleted together after merging.

As shown in Figure 5, step C further includes:

Step C1: Open the sequence file and read each sequence one by one;

Step C2: read the sequence (including the match of its reverse sequence) of each matching sequence of the variable-length kmer set, and count each variable-length kmer;

Step C3: traverse each bidirectional edge e in the De Bruijn graph;

Step C4: Count the maximum length of the sequencing sequences of all input files, and assign it to Lmax;

Step C5: If the length of the two-way side is greater than the length of Lmax-kmer, return to step C3, otherwise execute step C6;

Step C6: The starting vertex of the bidirectional edge e is u, the starting direction is Direc_u, the ending vertex is v, and the ending direction is Direc_v;

Step C7: Take all the edges that enter the vertex u and end in the direction of Direc_u, take the last k+1 characters of the last two-way edge e_i and put them into the string array m, and take all the two-way edges e_j that leave the vertex v and start in the direction of Direc_v The last character is put into the string array n;

Step C8: Combine the left and right sides of the two-way sides with the strings in the m and n arrays to construct a variable-length kmer array, wherein each element in the variable-length kmer array is variable-length kmer=m+e+n. Query these variable-length kmers in the counted variable-length kmers, and take out the one with the largest count value, and merge the corresponding three sides e_i, e, and e_j, and then delete e_i, e, and e_j.

When there are still two-way edges to be processed, return to step C3; when there are still sequences to be processed, return to step C1; when there are no two-way edges to be processed and no sequences to be processed, end.

Compared with prior art, the beneficial effect of the present invention is:

(1) Construct a variable-length kmer combination only for the existing bifurcation edges on the De Bruijn graph, and query the number of occurrences in the input sequence, and then select the merge of two-way edges according to the number; while IDBA iterates all kmer length, it is necessary to construct all possible kmers of each kmer length, and then shrink the De Bruijn graph, which method will lead to greater memory consumption and calculation time consumption;

(2) Select the optimal merging of bifurcated edges to minimize the possibility of erroneous merging of bidirectional edges;

(3) The length of contigs can be significantly increased, and the loss of contig quality can also be minimized; compared with other existing methods that increase the contig length and sacrifice contig quality, the present invention has a certain degree of control and improvement.

Those of ordinary skill in the art can understand that all or part of the steps in the various methods of the embodiments can be completed by instructing related hardware through a program. The program can be stored in a computer-readable storage medium, and the storage medium can include: only Read memory (ROM, Read Only Memory), random access memory (RAM, Random AccessMemory), disk or CD, etc.

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention should be included in the protection of the present invention. within range.

Claims (7)

1. a kind of two-way side extended method that elongated kmer based on the two-way De Bruijns of multistep is inquired about, it is characterised in that bag Include：

Step A：Read sequencing data source file, the two-way De Bruijns of construction multistep；

Step B：Count to intersecting two-way side in the two-way De Bruijns of the multistep；

Step C：The two-way side extension inquired about based on elongated kmer in the two-way De Bruijns of the multistep；

In step B, the two-way side of all the right and left summits bifurcated in De Bruijns two-way to the multistep and institute State it is two-way while bifurcated while the be possible to situation that merges, be assembled into a longer kmer, and inquire about correspondingly longer The occurrence number of kmer, the size with the occurrence number carry out selection and the union operation of bifurcated side trend as weight；It is described In step C, after all possible bifurcated on the two-way side to each bifurcated is converted to longer kmer, according to it is described go out occurrence Several height, will appear from number of times highest diverging paths as final bifurcated Merge Scenarios, and by corresponding diverging paths Two-way side delete together after merging.

2. the method for claim 1, it is characterised in that step B is further included：

Step B1：Each two-way side e in traversal De Bruijns；

Step B2：The maximum length of the sequencing sequence of all of input file is counted, and is assigned to Lmax；

Step B3：If the length on two-way side is more than the length of Lmax-kmer, execution step B1, otherwise execution step B4 are returned；

Step B4：The initial vertex of two-way side e is u, and prime direction is Direc_u, and termination summit is v, terminates direction and is Direc_v；

Step B5：All entrance summit u, terminate direction for Direc_u while take its it is two-way while last k+1 character be put into word Symbol string array m, it is all to leave the two-way side that vertex v, prime direction are Direc_v and take its last character and be put into character string number Group n；

Step B6：By two-way side or so respectively with m and n arrays in character string all combinations of construction, construct elongated kmer numbers Group, wherein each element in elongated kmer arrays is elongated kmer=m+e+n.

3. the method for claim 1, it is characterised in that step C is further included：

Step C1：Sequencing sequence file is opened, every sequence is read one by one；

Step C2：The sequence that elongated kmer collection charge-coupled each matching is read in, including the matching to its reverse sequence, and to each change Long kmer is counted；

Step C3：Each two-way side e in traversal De Bruijns；

Step C4：The maximum length of the sequencing sequence of all of input file is counted, and is assigned to Lmax；

Step C5：If the length on two-way side is more than the length of Lmax-kmer, execution step C3, otherwise execution step C6 are returned；

Step C6：The initial vertex of two-way side e is u, and prime direction is Direc_u, and termination summit is v, terminates direction and is Direc_v；

Step C7：All entrance summit u, terminate direction for Direc_u while take its it is last two-way while e_i last k+1 word Symbol is put into character string dimension m, all to leave the two-way side e_j that vertex v, prime direction are Direc_v and take its last character It is put into character string dimension n；

Step C8：By two-way side or so respectively with m and n arrays in character string all combinations of construction, construct elongated kmer numbers Group, wherein each element in elongated kmer arrays is elongated kmer=m+e+n, inquires about these changes in the elongated kmer for counting Long kmer, and maximum one of count value is taken out, and by corresponding e_i, tri- sides of e, e_j merge, and then delete e_i, E, e_j.

4. method as claimed in claim 3, it is characterised in that elongated kmer collection charge-coupled each matching is read in step C2 Sequence include the matching to reverse sequence.

5. the method for claim 1, it is characterised in that step A is further included：

Compression storing step, specially

A11, one sequence s of reading；

A12, sequence s sliding window is cut into into multiple fragments t；

A13, to each fragment t, encoded using nucleic acid coding table, and be expressed as the integer a of 64；

A14, fragment t is inverted, the fragment complementation for inverting is processed using symmetrical complement table, obtained complementary fragment v, and again Complementary fragment is encoded by the nucleic acid coding table in secondary use step A13, and is expressed as the integer b of 64；

The maximum number of A15, round numbers a and integer b, as the conventional number of the k molecules of fragment t and complementary fragment v；

A16, repeat step A11-A15, until all sequences are completed；

With De Bruijn constitution steps, specially

A21, one sequence s of reading；

A22, sequence s sliding window is cut into multiple fragments t, choose fragment t its conventional number be cur and before marking which, The conventional number of fragment afterwards is respectively pre, lat；

If the coding of A23, t is encoded less than its complementary fragment, pre, the value of lat are exchanged；

A24, cur forward position mapping table corresponding bit positions 1 come represent point to pre side；

A25, cur reverse position mapping table corresponding bit positions 1 come represent point to lat side；

A26, repeat step A22-A25, process other fragments t of sequence s, until completing whole fragments t of sequence s, perform step Rapid A27；

A27, one new sequence s of reading, repeat step A22-A26；Until having processed all of sequence, execution step A28；

A28, the construction for completing two-way multistep De Bruijn.

6. method as claimed in claim 5, it is characterised in that the sliding window in step A12, A22 is that length is k's Sliding window, wherein 0<k<32 and k is odd number；

Nucleic acid coding table in step A13 is { A:00,C:01,G:10,T:11}；

Symmetrical complement table in step A14 is { A->T,C->G,G->C,T->A}；

Step A14 specifically, the character string of fragment t is inverted, using in the character string that symmetrical complement table will be inverted Each character is changed into its complementary character, obtains character string v of complementary character, and reuses the nucleic acid coding table in step A13 Character string v is encoded, and is expressed as the integer b of 64.

7. method as claimed in claim 5, it is characterised in that in step A22, if before or after fragment t does not have Fragment, then be assigned to sky or NULL to pre or lat values；

In step A24, forward position mapping table is { A:0, C:1, G:2, T:3 }, the last character of position enquiring character for pre Symbol；

In step A25, reverse position mapping table is { A:4, C:5, G:6, T:7 }, first character of the position enquiring character for lat Complementary character.