A gap-filling method based on third-generation sequencing data

Abstract

An embodiment of the present invention discloses a gap-filling method based on third-generation sequencing data, comprising the following steps: obtaining an assembled sequence, dividing the assembled sequence into multiple overlapping groups according to gaps; aligning the multiple overlapping groups with each other, and removing excessively long overlapping sequences at the ends of the multiple overlapping groups; aligning reads with the overlapping group sequences, retaining reads that successfully align with the ends of the overlapping group sequences and reads that fail to align with the overlapping group sequences, with the retained reads being referred to as UTreads; aligning the UTreads with themselves, finding dovetail overlap relationships between reads, and deleting dovetail overlap relationships and reads with high depth; constructing a directed graph with the ends of the overlapping groups and both ends of the UTreads as nodes based on the alignments between the ends of the multiple overlapping groups, the alignments between the reads and the ends of the overlapping group sequences, and the dovetail overlap relationships obtained by aligning the UTreads with themselves; finding an optimal path between the ends of the overlapping groups on both sides of the gap based on the directed graph with the ends of the overlapping groups and both ends of the reads as nodes; and generating a sequence to fill the gap based on the optimal path.

Description

Technical Field

[0001] The present invention relates to the field of gene sequence gap filling technology, and more particularly to a gap filling method based on third-generation sequencing data. Background Art

[0002] The emergence of single-molecule sequencing (SMS) technology, along with a wide range of assembly software based on SMS-generated long reads, has significantly improved the continuity of genome assemblies for complex organisms, such as those found in eukaryotes. SMS technologies primarily include sequencing technologies from PacBio (Pacific Biosciences) and ONT (Oxford Nanopore Technology). Commonly used SMS long-read assembly software include Canu, Falcon, wtdbg2, NECAT / MECAT, Flye, Shasta, and NextDenovo. Further improvements in sequencing technology, such as the generation of ultra-long ONT sequencing reads, have enabled the automated assembly of human centromere regions. The high-fidelity reads generated by CCS (Circular Consensus Sequencing) sequencing have greatly improved the accuracy of SMS long reads and reduced the impact of sequencing errors on the assembly of highly similar repetitive regions within the genome. BioNano and Hi-C sequencing technologies assemble contigs into scaffolds or map them to chromosomes. If gaps within scaffolds can be filled using software and sequencing reads, assembly continuity and contig N50 can be further improved. The combination of these sequencing technologies and software enables the assembly of more complex genomes, such as those of eukaryotes, to produce genome reference sequences with improved continuity and more accurate repeat regions. Existing gap-filling software based on third-generation sequencing data includes TGS-GapCloser, LR_Gapcloser, Cobbler, GMcloser, FGAP, and PBJelly.

[0003] Existing gap-filling software based on third-generation sequencing data include TGS-GapCloser, LR_Gapcloser, Cobbler, GMcloser, FGAP, PBJelly, etc.

[0004] TGS-GapCloser software is a software that connects the two sides of the gap by aligning the long reads of TGS or another assembled sequence to the sequences on both sides of the gap. The input data of TGS-GapCloser are the long reads of TGS, the short reads of NGS sequencing, and the assembled scaffold sequence of the hole to be filled. The operation process is: detect the position of the gap in the scaffold sequence; use minimap2 to align the long reads to the sequences on both sides of the gap, extract the sequences aligned with the 2kb on both sides of the gap and the gap region as candidate fragments, and select the 10 candidate fragments with the highest QS scores for each gap; align the short reads to the candidate fragments, and use Pilon to correct the candidate fragments; use Racon to correct the candidate fragments based on the alignment of the long reads; re-align the candidate fragments to the sequences on both sides of the gap to calculate the QS score, and select the candidate fragment with the highest score to fill the gap sequence. See the detailed flowchart. Figure 2 However, this technology only uses single reads to fill the gap region, and cannot connect multiple reads to fill the gap in long gap regions.

[0005] LR_Gapcloser software generates 300bp tags based on long TGS reads and aligns the tags to the sequences flanking the gap to fill gaps. TGS-GapCloser takes long TGS reads and the assembled scaffold sequence to be filled. The program operates as follows: Long reads are split into 300bp non-overlapping tags. Tags are aligned to the assembled scaffold sequence using BWA-MEM, retaining alignments with greater than 80% coverage. For tags that uniquely align to a specific scaffold position, the tag with the most tags aligned is retained. The median tag is then selected, and the tag whose left and right sides are aligned in the same order in the scaffold is retained, using this tag as a reference. For tags that align to multiple scaffold positions, alignment positions are selected based on alignment direction consistency with adjacent uniquely aligned tags and alignment order in the read that is consistent with the scaffold. If the tag in a read is aligned to both sides of the gap at the same time, and the ratio of the tag coordinate difference in the read to the coordinate difference in the scaffold meets a specific condition, then the tag pair is retained. Each gap is filled based on the read sequence with the most tag pairs. If the tag of a read is aligned to only one side of the gap, one side of the gap is extended based on the longest read covering the gap area. The entire gap filling process can be repeated multiple times to fill as many gaps as possible. See the detailed flowchart for details. Figure 3However, this technology only uses a single read to fill the gap area in each gap filling process, and cannot use multiple reads to connect and fill the gap in long gap areas. Summary of the Invention

[0006] To this end, an embodiment of the present invention provides a gap-filling method based on third-generation sequencing data to solve the problem that the existing technology cannot fill gaps in long gap regions.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for filling holes in third-generation sequencing data, comprising the following steps:

[0008] Obtaining an assembled sequence, and dividing the assembled sequence into multiple contigs according to gaps;

[0009] Aligning the plurality of contigs to each other, and removing excessively long overlapping sequences at the ends of the plurality of contigs;

[0010] Align the reads with the contig sequence, retain the reads that successfully align with the ends of the contig sequence and the reads that fail to align with the contig sequence, and the retained reads are called UT reads;

[0011] The UT reads were aligned with each other to find the dovetail overlap relationships between reads, and high-depth dovetail overlap relationships and reads were deleted;

[0012] Constructing a directed graph with the contig ends and both ends of the UT reads as nodes based on the dovetail overlap relationships obtained by alignment between the ends of the contigs, alignment between the reads and the ends of the contig sequences, and alignment of the UT reads themselves;

[0013] Finding the optimal path between the ends of the contigs on both sides of the gap based on a directed graph with the ends of the contigs and the ends of the reads as nodes;

[0014] A sequence gap is repaired according to the optimal path.

[0015] Furthermore, the method of obtaining the assembled sequence and dividing the assembled sequence into multiple overlapping groups according to the gap specifically includes: comparing the overlapping groups according to parameters, finding the parts with dove-tail overlapping relationships between the ends of the overlapping groups, counting the lengths of the ends of the overlapping group sequences covered, and deleting the overlapping group when the entire length of a certain overlapping group is covered; if the length of the end of a certain overlapping group covered is greater than a preset threshold, deleting the end of the overlapping group that exceeds the threshold coverage area.

[0016] Furthermore, the reads are aligned with the contig sequence, and the reads that are successfully aligned with the ends of the contig sequence and the reads that are not successfully aligned with the contig sequence are retained. The retained reads are called UT reads, and specifically include:

[0017] The reads are aligned to the contig sequences, reads aligned to the middle of the contig are removed, and reads aligned to the ends of the contigs and not aligned to the contigs are retained. The retained reads aligned to the ends of the contigs and not aligned to the contigs are UT reads.

[0018] Furthermore, the UT reads are aligned with each other to find the dovetail overlap relationships between reads, and high-depth dovetail overlap relationships and reads are deleted, specifically including:

[0019] UT reads are aligned to themselves to find alignments with dove-tail overlaps, and reads and alignments in high-depth regions are filtered out.

[0020] Furthermore, the dove-tail overlapping relationship obtained by aligning the ends of the multiple overlapping groups, aligning the reads with the ends of the overlapping group sequences, and aligning the UT reads themselves is used to construct a directed graph with the ends of the overlapping groups and the two ends of the UT reads as nodes, specifically including: constructing a directed graph with the ends of the overlapping groups and the two ends of the reads as nodes according to the dove-tail overlapping relationship obtained by aligning the ends of the multiple overlapping groups, aligning the reads with the ends of the overlapping group sequences, and aligning the UT reads themselves, deleting useless nodes and high-depth nodes that will not be accessed when finding a path, retaining a certain number of outgoing edges and incoming edges for the end nodes and reads nodes of the overlapping groups, respectively, and deleting the corresponding incoming edges of a portion of the nodes that form a loop without affecting the connectivity between all the overlapping group ends of the graph.

[0021] Furthermore, the step of finding the optimal path between the ends of the contigs on both sides of the gap based on a directed graph with the ends of the contigs and the ends of the reads as nodes specifically includes:

[0022] Find the path between the ends of the overlapping groups on both sides of the gap. The specific path finding method includes one of Dijkstra algorithm path finding, traversal path finding and step path finding.

[0023] Each of the Dijkstra algorithm pathfinding, traversal pathfinding, and step pathfinding methods screens edges in a directed graph based on the gradient of similarity and alignment length scores, and finds multiple paths.

[0024] Furthermore, the Dijkstra algorithm pathfinding will extract a subgraph starting from a contig end node on one side of the gap and terminating at the corresponding contig end node on the other side of the gap, and use the Dijkstra algorithm to find multiple optimal paths from the aforementioned starting node to all predecessor nodes of the aforementioned terminating node, and the optimal paths found are recorded according to the similarity and comparison length score intervals; each time the candidate extension node is calculated in the Dijkstra algorithm, the minimum score of all edges in the path from the starting node to the current candidate extension node is calculated, and each time the extension is performed, the end node of the path with the largest minimum score is selected to extend forward, and the final path obtained is the path with the largest minimum score of all edges in the path from the starting node to a certain terminating node.

[0025] Furthermore, finding the optimal path between the ends of the contigs on both sides of the gap based on the directed graph between the ends of the contigs and the ends of the reads specifically includes:

[0026] Each path is segmented by path extension length, and a path score sum is calculated based on the path extension length segmentation. Based on the distribution of the path score sum within each length segment, a distribution peak that meets the conditions is selected, and the path with the highest score within the peak is selected.

[0027] A method for filling a hole in a sequence using an assembled sequence to complete another assembled sequence based on third-generation sequencing data, comprising the following steps:

[0028] Extract the sequences on both sides of the gap of the first assembled sequence and align them to the second assembled sequence;

[0029] Filtering the alignment of the first assembled sequence and the second assembled sequence and generating a backbone;

[0030] Filter the backbone according to the length of the unaligned regions at both ends of the sequence;

[0031] Extracting original alignment information based on the main chain to generate a path;

[0032] A hole-filled sequence is generated according to the path.

[0033] Embodiments of the present invention have the following advantages: They can fill larger gaps using long sequences formed by overlapping multiple reads. For gaps, sequences on either side of the gap and potential reads in the gap region are aligned to obtain overlaps. A directed graph is constructed to find paths between the sequences on either side of the gap, and the sequences along the paths are used as sequences spanning the gap region. Another assembled sequence can be used to fill gaps. If another assembled sequence contains sequences spanning certain gap regions, the sequence in the gap region from the other assembled sequence is used to fill the gaps. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for the embodiments or the description of the prior art. Obviously, the drawings described below are merely exemplary, and those skilled in the art can derive other implementation drawings based on the provided drawings without inventive effort.

[0035] The structures, proportions, sizes, etc. illustrated in this specification are intended only to complement the contents disclosed herein and to facilitate understanding and reading by persons familiar with the art. They are not intended to limit the conditions under which the present invention may be implemented and therefore have no substantive technical significance. Any structural modifications, changes in proportions, or adjustments in sizes, without affecting the efficacy and objectives of the present invention, shall still fall within the scope of the technical contents disclosed herein.

[0036] Figure 1 A schematic diagram of a method flow for a hole filling method based on third-generation sequencing data provided by an embodiment of the present invention;

[0037] Figure 2 A flowchart of the TGS-GapCloser software for a gap-filling method based on third-generation sequencing data provided by an embodiment of the present invention;

[0038] Figure 3 LR_Gapcloser software flow chart for a gap filling method based on third-generation sequencing data provided by an embodiment of the present invention;

[0039] Figure 4 A flowchart of a gap-filling method based on third-generation sequencing data provided by an embodiment of the present invention using read sequences to fill gaps in assembled sequences;

[0040] Figure 5 A flowchart of a hole-filling method based on third-generation sequencing data provided by an embodiment of the present invention, which uses another assembly sequence B to fill holes in a certain assembly sequence A. DETAILED DESCRIPTION

[0041] The following describes the implementation of the present invention using specific embodiments. Those skilled in the art will readily understand the other advantages and benefits of the present invention from the disclosure herein. Obviously, the embodiments described are only a portion of the present invention, not all of it. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without inventive effort are intended to fall within the scope of protection of the present invention.

[0042] For a better understanding, the following professional terms are explained and annotated:

[0043] Third-generation sequencing (TGS) is a long-range DNA sequencing technology. Representative technologies include Pacific Biotechnology's SMRT technology in the United States and Oxford Nanopore's nanopore technology in the United Kingdom. For example, nanopore sequencing involves direct detection of the electrical signal generated by each DNA molecule as it passes through a nanopore, without PCR amplification, and converting it into a base sequence for real-time, single-molecule de novo sequencing.

[0044] Next Generation Sequencing (NGS): A short-range DNA sequencing technology. Representative technologies include Illumina's HiSeq sequencing technology and MGI's CG sequencing technology. For example, HiSeq sequencing technology uses a sequencing-by-synthesis approach. In template-based PCR sequencing, each round of synthesis reaction adds a fluorescently labeled base to the end of the sequence. After scanning and reading the base, the protective base is removed and the next round of reaction continues until sequencing is complete.

[0045] Genome: refers to the sum of all the genetic material DNA of an organism.

[0046] Sequencing sequence (read): A DNA molecular sequence information consisting of four bases, AGTC, obtained through DNA sequencing technology.

[0047] DNA sequence alignment: Two DNA sequences are aligned based on the order of the AGTC bases they contain, with matches, mismatches, gaps, etc. Two DNA sequences are considered aligned only if there are enough matches.

[0048] Overlap: If sequencing sequence A and sequencing sequence B can be aligned and there is an overlapping region, then sequencing sequence A and sequencing sequence B are said to have an overlapping relationship.

[0049] Dovetail overlap: If the overlap between sequencing sequence A and sequencing sequence B is located at the ends of the two sequencing sequences, and the length of the uncovered region at the ends is less than the threshold, the overlap between sequencing sequence A and sequencing sequence B is called a dovetail overlap.

[0050] Directed graph: A graph consists of a nonempty finite set of vertices and a finite set of edges. Let V and E represent the vertex set and edge set of graph G, respectively. Then, graph G can be represented as G = (V, E). If each edge of a graph has an arrow indicating its direction, it is a directed graph.

[0051] A directed graph (String Graph) is constructed based on the dovetail overlap relationship of sequenced data. A vertex, also called a node, represents an end of a sequenced or assembled data. Edges represent dovetail overlap relationships between sequenced and / or assembled data.

[0052] Contig: A continuous sequence without gaps obtained by assembly.

[0053] N50: An indicator of sequence continuity. Sort the sequences from longest to shortest and add up their lengths. When the sum reaches a certain sequence, and the accumulated length reaches 50% of the total sequence length, the length of that sequence is the N50 length of the group.

[0054] Example 1: A method for filling holes based on third-generation sequencing data, such as Figure 1 As shown, the following steps are included:

[0055] Obtain the assembled sequence and divide the assembled sequence into multiple contigs according to the gap;

[0056] Align multiple contigs with each other and remove the excessively long overlapping sequences at the ends of multiple contigs;

[0057] The reads were aligned with the contig sequences, and the reads that were successfully aligned with the ends of the contig sequences and the reads that were not successfully aligned with the contig sequences were retained. The retained reads were called UT reads;

[0058] Align UT reads with themselves to find dovetail overlaps between reads and delete high-depth dovetail overlaps and reads.

[0059] Based on the dove-tail overlap relationships obtained by alignment between multiple contig ends, alignment between reads and contig sequence ends, and alignment of UT reads themselves, a directed graph (digraph) with contig ends and both ends of UT reads as nodes was constructed;

[0060] Based on a directed graph with the ends of the contigs and the ends of the reads as nodes, find the optimal path between the ends of the contigs on both sides of the gap;

[0061] Generate a sequence to fill the gap based on the optimal path.

[0062] Specifically, refer to Figure 1 and Figure 4As shown, (1) the contigs are divided according to the gap, and the overlap at the ends of the contigs is removed. The scaffold sequence is divided into contigs according to the gap, and the contig name information on both sides of the gap is recorded. The contigs are aligned according to the parameters (software and parameters: minimap2-x asm5-DP-f 100-g 100000-r 10000-n 10--dual=no) to find the dovetail overlap between the ends of the contigs. The length of the ends of the contig sequences covered is counted. If the entire length of the contig is covered, the contig is deleted. If the length of a certain end covered is >25kb, the end of the contig with the extra 25kb coverage is deleted. According to the length of the deleted end, the dovetail alignment relationship between the contigs is adjusted.

[0063] (2) Generate UT reads and align them to themselves. Align the reads to the genome, remove reads aligned to the middle of the contig, and retain reads aligned to the ends of the contig and those not aligned to the contig (UT reads). Alignment software and parameters for hifi reads: minimap2-x asm5-a (the new version of minimap2 can also use the -xmap-hifi-a parameter). Alignment software and parameters for ontreads: minimap2-x map-ont-a. Based on the alignment of reads aligned to the ends of the contig, generate the dovetail overlap between the end of the contig and one end of the read. Align the UT reads to themselves, find the dovetail alignment, and filter out reads and alignments in high-depth areas (average depth > 1000). Alignment software and parameters for hifi reads: minimap2-x asm5-XHk 51-w 51-n 10-f 1000 (minimap2 software has been modified to support k being set to any value in the range [1,128]). Ont reads self-alignment software and parameters: minimap2 -x ava-ont-k17 -w 17 -n 10 -f 1000.

[0064] (3) Construct a directed graph of the contig ends and the two ends of the reads. Based on the alignment between contigs, the alignment between the contig ends and one end of the reads, and the alignment between the dovetails of the reads, construct a directed graph (string graph) with the contig ends and the two ends of the reads as nodes. Delete useless nodes that will not be accessed by the path finding. Delete high-depth nodes (nodes with an out-degree greater than a threshold, where the threshold is 10 times the median out-degree of all nodes). For the contig end nodes and read nodes, retain a certain number of outgoing edges and incoming edges respectively. Without affecting the connectivity between all contig ends of the graph, delete the corresponding incoming edges of some nodes that form a loop (this step is only required for traversal path finding and Dijkstra algorithm path finding).

[0065] (4) Find the path between the ends of the contigs on both sides of the gap. There are three methods, as follows.

[0066] Traversal path search: Extract the subgraph starting from the end node of a contig and ending at the end node of a contig, breaking all loops in the subgraph. Traverse all nodes using a breadth-first search algorithm, record the best path between intervals based on similarity and dt_score, and calculate the path score avg_score, etc.

[0067] Dijkstra's algorithm for finding paths (default method): Extract a subgraph starting from the end node of a contig on one side of the gap and terminating at the end node of the corresponding contig on the other side of the gap. Use the Dijkstra algorithm to find the path with the highest edge score each time it is extended. The best path is recorded between partitions based on similarity and dt_score. The detailed process is as follows: Nodes are divided into two sets, U and V, representing nodes with and without distances calculated, respectively. Initially, U contains only the start node S, and V contains all remaining nodes. Each time, from edges in U pointing to nodes in V, the edge with the highest score is selected and extended. The terminating node of this edge is removed from V and added to U. Repeat this process until paths from S to all other nodes are found. Each time, the minimum score of all edges in the path is calculated, and the edge with the highest minimum score is selected for extension. The resulting path is the path with the highest minimum score among all edges in the path from the start point to each node. The scoring function is score = (dt_score + 300) / 300. Where dt_score = align_len / 5 - (left_overhang + right_overhang) / 2. align_len is the overlap length of the edge alignment, and left_overhang and right_overhang are the lengths of the unaligned regions on one side of the dovetail overlap alignment between the two reads. Calculate the paths and path scores between all predecessor nodes from the start node to the end node. Merge multiple paths into a single path. The avg_score of the path is the average of the scoring function values, and the number of paths is the number of merged paths.

[0068] Find paths in a stepwise manner: Filter edges by setting a minimum threshold based on similarity and dt_score partitioning. Starting from a contig end node, select nodes based on the maximum dt_score (best_score), maximum extension length (max_extending), and randomization. Extend the path forward until the next contig end node is found. Record the found path based on similarity and dt_score partitioning. Merge multiple paths into a single path. The avg_score of the path is the matches percentage value of each edge in the path - 80. The prop_diff_sum value is then summed up. The prop_diff_sum values of the multiple paths are then averaged to obtain the avg_score. The number of paths is the number of paths.

[0069] (5) Select a path. For the path found in the previous step, calculate the length segment index number based on the path extension length. Lengths in the range [0, 5000) are segmented into 500bp segments with an index range of [0, 9]; lengths in the range [5000, 100000) are segmented into 2kb segments with an index range of [10, 57]; lengths in the range [100000, 500000) are segmented into 10kb segments with an index range of [58, 107]; lengths in the range [500000, 5000000) are segmented into 100kb segments with an index range of [108, 152]; lengths in the range [5000000, 50000000) are segmented into 1Mb segments with an index range of [153, 197]; negative lengths are segmented into 500bp segments with an index range of [-inf, 0].

[0070] The sum of the path scores for each length segment is calculated. The distribution is calculated based on the sum of the path scores within each length segment. A peak is defined as a score greater than or equal to the scores within the segments on both sides, while a valley is defined as a score less than the scores within the segments on both sides, or less than the scores within at least one segment on either side. If the score within a valley is less than 50% of the score within the next peak, the valley is split into multiple peaks. Within each peak, segments with scores greater than or equal to 0.1 times the highest score are retained, and the sum of the scores of the retained segments within the peak is calculated. The peak with the highest sum of scores within the peak is selected. If the score of the next longer peak is greater than or equal to 50% of the score of the previous peak, the next peak is selected, and this process is repeated until the last peak is selected. For each selected peak, the path with the largest sum of ovl_scores within the peak is selected as the final path.

[0071] The ovl_score of an edge = (query_align_len + target_align_len) * identity - left_overhang - right_overhang. Path score = path type weight * avg_score * number of paths * weights of the minimum path similarity and minimum dt_score. Path type weights: Dijkstra paths have a weight of 3, best_score paths for stepwise paths have a weight of 3, max_extending paths have a weight of 1, and randomizing paths have a weight of 0.5. For avg_score and the number of paths, see the description of the path finding method above.

[0072] The weights for the minimum path similarity and minimum dt_score are calculated based on the ranking of the path's minimum similarity among all path minimum similarities and the ranking of the path's minimum dt_score among all path minimum dt_scores. The maximum and minimum values of all path minimum similarities are calculated, divided into 100 windows, and the sum of the number of paths in each window is counted. The windows are sorted from largest to smallest by index, and a mid-similarity threshold is set. The total number of paths greater than the mid-similarity threshold, mid_num, is calculated. The windows are traversed from largest to smallest, and the total number of paths is added to total_num. If the similarity of the left edge of the window is greater than or equal to the mid-similarity threshold, the weight of the similarity within that window is pow(mid_num, 0.2) / pow(total_num, 0.2); otherwise, the weight is pow((mid_num + 1) / total_num, 0.9). The weight for the minimum dt_score is calculated similarly.

[0073] (6) Generate the gap-filled sequence. Based on the path selected in the previous step, generate the sequence coordinates of each read to be intercepted in the path, merge the intercepted sequences of each read to obtain the sequence of the gap region, replace the original contig end sequence, and connect the corresponding contigs.

[0074] The following example uses a test bacterium with multiple gaps to obtain the scaffold sequence of the test bacterium with multiple gaps, as well as the Pacbio sequencing reads sequence used to fill the gaps.

[0075] The bacterial scaffold sequence was divided into contig sequences according to the gap, and the contig name information on both sides of the gap was recorded. The contigs were aligned according to the aforementioned parameters to find the dovetail overlap between the contig ends. The length of the covered contig sequence ends was counted. If the entire contig length was covered, the contig was deleted. If the covered length of a certain end was >25kb, the contig end with the excess coverage of 25kb was deleted. The dovetail alignment relationship between contigs was adjusted according to the length of the deleted end. Contig sequences <50kb were used as reads, and contig sequences >=50kb were used as contig sequences for gap filling.

[0076] Reads were aligned to the genome using the aforementioned hifi reads parameters, with reads aligned to the middle of the contig removed. Reads aligned to the ends of the contig and reads not aligned to the contig (UT reads) were retained. Based on the alignment of reads to the ends of the contig, dovetail overlaps were generated between the end of the contig and one end of the read.

[0077] The UT reads were aligned to themselves according to the aforementioned hifi reads parameters to find dovetail alignments, and reads and alignments in high-depth regions (average depth > 1000) were filtered out.

[0078] Based on the alignment between contigs, the end of a contig and one end of a read, and the dovetail alignment between reads, a string graph is constructed with the ends of the contig and both ends of the read as nodes. Unused nodes that will not be accessed by the path search are deleted. High-depth nodes (nodes with an out-degree greater than a threshold, where the threshold is 10 times the median out-degree of all nodes) are deleted. For contig end nodes and read nodes, a certain number of outgoing and incoming edges are retained (the 40 highest-scoring outgoing and incoming edges are retained for contig end nodes, and the 20 highest-scoring outgoing and incoming edges are retained for read nodes). Some incoming edges corresponding to nodes that form a loop are deleted.

[0079] A path is found using the Dijkstra algorithm. Each time, the minimum score of all edges in the path is calculated. Each time an edge is extended, the path with the highest minimum score is selected for forward extension. Paths and path scores are calculated for all predecessors from the start node to the end node. Multiple paths are merged into a single path. The scoring function is described above. For each found path, the length segment index is calculated based on the path's extended length. The sum of the path scores for each length segment is calculated. Based on the sum of the path scores within each length segment, the peaks and valleys of the distribution are determined. The distribution is divided into multiple peaks based on the valleys. Within each peak, segments with a score greater than or equal to 0.1 times the highest score are retained. The sum of the scores of the retained segments within the peak is calculated. The peak with the highest sum of scores within the peak is selected. If the next longer peak has a score greater than or equal to 50% of the previous peak's score, the next peak is selected. This process is repeated until the final peak is selected. For each selected peak, the path with the highest sum of ovl_scores for the edges within the peak is selected as the final path. The scoring function is described above. According to the selected path, the sequence coordinates of each read to be intercepted in the path are generated, the sequence of each read interception is merged to obtain the sequence of the gap region, the original contig end sequence is replaced and the corresponding contigs are connected.

[0080] Large gaps can be filled using long sequences formed by overlapping multiple reads. For gaps, the overlap is obtained by aligning the sequences on either side of the gap and the potential reads in the gap. A directed graph is constructed to find the path between the sequences on either side of the gap, and the sequence along the path is used as the sequence spanning the gap.

[0081] Example 2: A method for filling a hole in a sequence using a certain assembly sequence to complete another assembly sequence based on third-generation sequencing data, referring to Figure 5 As shown, the following steps are included:

[0082] Extract the sequences on both sides of the gap of the first assembled sequence and align them to the second assembled sequence;

[0083] The alignment of the first assembled sequence and the second assembled sequence is filtered and a backbone is generated;

[0084] Filter the backbone based on the length of the unaligned regions at both ends of the sequence;

[0085] Extract the original alignment information based on the main chain and generate a path;

[0086] Generate a hole-filled sequence based on the path.

[0087] Specifically, (1) extract the sequences on both sides of the gap of the assembled scaffold sequence A and align them to the assembled sequence B. The software and parameters for aligning the sequences on both sides of the gap to the assembled sequence B are: minimap2-x asm5-I 16G-f 100.

[0088] (2) Filter the alignments and generate the backbone. Remove alignments with alignment lengths < 100k (1 / 10 of the default length parameter 1Mb). Sort the starting coordinates of the query sequence (sequences on both sides of the gap) and the target sequence (sequence in the assembled sequence B) by the aligned fragments. If the gap or overlap between adjacent fragment pairs is less than a threshold, connect the adjacent aligned fragment pairs to form the backbone. The gap length threshold in A is 1Mb (default length parameter 1Mb), the gap length threshold in B is 2Mb (twice the default length parameter 1Mb, which controls the maximum length of the gap-filling sequence in B), and the overlap length threshold is 500kb (1 / 2 of the default length parameter 1Mb). Remove backbones with lengths less than the threshold 500kb (1 / 2 of the default length parameter 1Mb).

[0089] (3) Filter the backbone based on the termini. If the unaligned regions on either side of the query sequence are less than 250 kb (1 / 4 of the default length parameter of 1 Mb), the backbone is retained.

[0090] (4) Extract the original alignment information based on the main chain and generate a path.

[0091] (5) Generate the gap-filled sequence. Based on the path information generated in the previous step, intercept the sequence of the corresponding coordinate region of each contig, merge the intercepted sequences to obtain the sequence of the gap region, replace the original contig end sequence, and connect the corresponding contigs.

[0092] Taking a test bacterium with multiple gaps as an example, a scaffold sequence A of the test bacterium with multiple gaps and a sequence B of another assembled version of the bacterium used to fill the gaps are obtained.

[0093] Extract the sequences on both sides of the gap in sequence A and align them to the assembled sequence B according to the aforementioned parameters. The default length parameter is set to 50kb. Remove alignments with alignment lengths <5k. Sort the starting coordinates of the query sequence (sequences on both sides of the gap) and the target sequence (sequence in assembled sequence B) according to the aligned fragments. If the gap or overlap between adjacent fragment pairs is less than the threshold, connect the adjacent aligned fragment pairs to form a main chain. The gap length threshold in A is 50kb, the gap length threshold in B is 100kb, and the overlap length threshold is 25kb. Remove the main chain with a length less than the threshold of 25kb. If the unaligned regions at both ends of the query sequence are <12.5kb, retain the main chain. Extract the original alignment information based on the main chain and generate a path. Based on the generated path information, intercept the sequence of the corresponding coordinate region of each contig, merge the intercepted sequences to obtain the sequence of the gap region, replace the original contig end sequence, and connect the corresponding contigs.

[0094] Another assembled sequence can be used to fill in gaps. If another assembled sequence has sequences that span the gaps in some gaps, the sequence in the gaps can be used to fill in the gaps.

[0095] Although the present invention has been described in detail above using general descriptions and specific embodiments, it will be apparent to those skilled in the art that modifications and improvements may be made thereto. Therefore, such modifications and improvements, which do not depart from the spirit of the present invention, are intended to fall within the scope of protection claimed herein.

Claims

1. A gap filling method based on third-generation sequencing data, characterized by: The following steps are involved: Obtaining an assembled sequence, and dividing the assembled sequence into a plurality of contig sequences according to gaps; Comparing the plurality of contig sequences with each other, and removing the overlapping sequences at the ends of the plurality of contig sequences whose lengths are greater than a preset length threshold; The reads are aligned with the contig sequence after removing the overlapping sequences at the ends of the multiple sequences with a length greater than a preset length threshold, and the reads that are successfully aligned with the ends of the contig sequence and the reads that are not successfully aligned with the contig sequence are retained. The retained reads are called UT reads; The UT reads are aligned with each other to find the dovetail overlap relationships between reads, and high-depth dovetail overlap relationships and reads are deleted, where the high depth is defined as the out-degree greater than a threshold; Constructing a directed graph with the ends of the contig sequences and the UT reads as nodes based on the dovetail overlap relationships obtained by alignment between the ends of the contig sequences, alignment between the reads and the ends of the contig sequences, and alignment of the UT reads themselves; Based on a directed graph with the contig sequence ends and UT reads as nodes, find the optimal path between the contig sequence ends on both sides of the gap; A sequence gap is repaired according to the optimal path.

2. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The plurality of overlapping group sequences are compared with each other, and overlapping sequences at the ends of the plurality of overlapping group sequences whose lengths are greater than a preset length threshold are removed, specifically comprising: comparing the overlapping group sequences according to parameters, finding the portions with dove-tail overlapping relationships between the ends of the overlapping group sequences, counting the covered lengths of the ends of the overlapping group sequences, and deleting the overlapping group sequence when the entire length of a certain overlapping group sequence is covered; if the covered length of the end of a certain overlapping group sequence is greater than a preset threshold, deleting the end of the overlapping group sequence that exceeds the threshold coverage area.

3. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The reads are aligned with the contig sequence after removing overlapping sequences with a length greater than a preset length threshold at the ends of multiple sequences, and the reads that are successfully aligned with the ends of the contig sequence and the reads that are not successfully aligned with the contig sequence are retained. The retained reads are called UT reads, specifically including: The reads are aligned with the contig sequence, reads aligned to the middle of the contig sequence are removed, and reads aligned to the end of the contig sequence and not aligned to the contig are retained. The retained reads aligned to the end of the contig sequence and not aligned to the contig are UT reads.

4. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The UTreads are aligned with each other to find dovetail overlap relationships between reads, and high-depth dovetail overlap relationships and reads are deleted, specifically including: UT reads are aligned with themselves, the parts with dove-tail overlap are found for alignment, and reads with overlap depth greater than the preset overlap depth threshold are filtered out.

5. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The method constructs a directed graph with the overlapping group sequence ends and the two ends of the UT reads as nodes based on the dove-tail overlapping relationship obtained by the alignment between the ends of the multiple overlapping group sequences, the alignment between the reads and the ends of the overlapping group sequences, and the alignment of the UT reads themselves, specifically including: constructing a directed graph with the overlapping group sequence ends and the two ends of the reads as nodes based on the dove-tail overlapping relationship obtained by the alignment between the ends of the multiple overlapping group sequences, the alignment between the reads and the ends of the overlapping group sequences, and the alignment of the UT reads themselves, deleting useless nodes and high-depth nodes that will not be accessed when finding a path, retaining a certain number of outgoing edges and incoming edges for the end nodes and reads nodes of the overlapping group, respectively, and deleting the corresponding incoming edges of a portion of the nodes that form a loop without affecting the connectivity between all overlapping group sequence ends of the graph.

6. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The directed graph based on the contig sequence ends and UT reads as nodes, finding the optimal path between the contig sequence ends on both sides of the gap, specifically includes: Find the optimal path between the ends of the contig sequences on both sides of the gap, and the path finding method includes one of Dijkstra algorithm path finding, traversal path finding, and stepping path finding; Each of the Dijkstra algorithm pathfinding, traversal pathfinding, and step pathfinding methods screens edges in a directed graph based on the gradient of similarity and alignment length scores, and finds multiple paths.

7. The method for filling holes in third-generation sequencing data according to claim 6, wherein: The Dijkstra algorithm pathfinding will extract a subgraph starting from a contig end node on one side of the gap and terminating at the corresponding contig end node on the other side of the gap, and use the Dijkstra algorithm to find multiple optimal paths starting from a contig end node on one side of the gap to the corresponding contig end node on the other side of the gap. The optimal paths found are recorded according to the similarity and comparison length score intervals; each time the candidate extension node is calculated in the Dijkstra algorithm, the minimum score of all edges in the path from the starting node to the current candidate extension node is calculated, and each time the extension is performed, the end node of the path with the largest minimum score is selected to extend forward, and the final path is the path with the largest minimum score of all edges in the path from the starting node to a certain end node.

8. The method for filling holes in third-generation sequencing data according to claim 1, wherein: The directed graph based on the contig sequence ends and UT reads as nodes, finding the optimal path between the contig sequence ends on both sides of the gap, specifically includes: Each path is segmented according to the path extension length, and the total path score is calculated based on the segmentation of each path according to the path extension length. According to the distribution of the total path score in each length segment, a distribution peak that meets the conditions is selected, and the path with the highest score in the peak is selected.

Images