Methods and systems for filling gaps in genome assembly

Abstract

The present application discloses a method and system for filling the gap of genome assembly. The present application filters out the conflicting alignment by processing the preliminary alignment file, and calculates a series of specific scores to evaluate the quality of the alignment, and obtains the sequence corresponding to the optimal gap filling based on the alignment quality. The present application can effectively fill the gap during genome assembly by simple operation without professional technical personnel, and is particularly suitable for telomere-to-telomere level genome assembly. Meanwhile, the present application can greatly reduce the gap filling time and avoid high workload.

Description

Technical Field

[0001] This invention relates to the field of genome assembly, and more particularly to methods for filling gaps in genome assembly, especially suitable for genome assembly processes from telomere to telomere level. Background Technology

[0002] In recent years, with the advancement of long-read sequencing technology, telomere-to-telomere (T2T) reference genomes for at least dozens of eukaryotes (such as humans and rice) have been established. This has not only promoted the structural and functional studies of repetitive genomic regions but also avoided errors in bioinformatics analysis, such as the detection of false positive variants and data contamination caused by incomplete reference genomes. However, since automated genome assembly processes cannot generate complete chromosome sequences, almost all of these projects include time-consuming manual gap-filling processes to supplement missing sequences from complex genomic regions (such as tandem and scattered repeats). In most T2T projects, automated assembly processes take only a few days, while manual processes can take months or even years. Furthermore, manual work requires high information skills, such as programming and Linux operation, which hinders the wider participation of biologists with limited information skills in T2T projects. In summary, due to the high workload and technical requirements, manual gap-filling has become a major obstacle to the T2T-level assembly of large eukaryotic genomes (such as wheat) and large-scale pan-genomes.

[0003] The information in the background section is merely intended to illustrate the general background of the invention and should not be construed as an admission or implication in any way that such information constitutes prior art known to those skilled in the art. Summary of the Invention

[0004] To address at least some of the technical problems in the prior art, the present invention provides a method and system for filling gaps in genome assembly. Specifically, the present invention includes the following.

[0005] A first aspect of the present invention provides a method for filling gaps in genome assembly, comprising the following steps:

[0006] (1) Provide chromosome-level assembly data and whole-genome third-generation sequencing data containing multiple read sequences generated from third-generation sequencing data;

[0007] (2) Generate a first alignment file and a second alignment file, wherein the first alignment file contains alignment data between reads and chromosomes, and the second alignment file contains alignment data between different reads;

[0008] (3) Filter out the read segments that do not have gap regions and retain the read segments with gap regions, wherein the gap regions include the gap region and the two sides of the gap;

[0009] (4) Process the first comparison file and the second comparison file to filter out conflicting comparison data and obtain valid comparison data;

[0010] (5) Evaluate the quality of each alignment in the valid alignment data obtained in step (4) and generate a reliability score;

[0011] (6) Optionally, a set of candidate reads for filling gaps is generated using a heuristic breadth-first search based on the reliability score.

[0012] (7) Input the chromosome-level assembly results, sequencing data, first and second alignment files and reliability score files into the personal terminal, and the personal terminal can visualize and select appropriate sequences to fill the gaps.

[0013] In some embodiments, the method for filling gaps in genome assembly according to the present invention further includes, in step (1), providing data for generating assembly units for chromosome-level assembly.

[0014] In some embodiments, the method for filling gaps in genome assembly according to the present invention includes the step of filtering out conflict alignments, which comprises:

[0015] In the fragment alignment set A between the reference sequence r and the query sequence q, the alignments (a i According to its corresponding segment on the reference sequence r (r i Sort the reference coordinates of the query sequence q in ascending order to obtain the corresponding segment (q) on the query sequence q. i ) sequence S q They are arranged in the same order, and the reference coordinates and S are obtained simultaneously. q q in the same order i Reference coordinates of sequence S c S is calculated using a binary search algorithm. c The longest increasing subsequence S c ', by collecting all segments from A that correspond to it on the query sequence q (q i The reference coordinates of ) are contained in S c The comparison within ' is used to obtain the optimal conflict-free subset A', thereby filtering out conflicting comparisons.

[0016] In some embodiments, according to the method for filling gaps in genome assembly according to the present invention, the reference coordinates include start coordinates and / or end coordinates, wherein the longest increasing subsequence S c The following condition must be met: for any two values ​​x and y in the sequence, if x comes before y, then x ≤ y.

[0017] In some embodiments, the method for filling gaps in genome assembly according to the present invention is wherein the quality of the alignment is assessed by calculating at least one of the following: conflict score CS, maximum alignment ratio PLA, total length of conflict-free alignments LNA, proportion of conflict-free alignments PNA, and number of conflict-free matching bases MBN.

[0018] In some embodiments, according to the method for filling gaps in genome assembly according to the present invention, the conflict score CS is calculated as follows:

[0019] A certain number of k-mers are randomly selected within the alignment region of the reference sequence, and their positions in the query sequence are calculated. The position of the k-mer in the query sequence is obtained based on the start and end coordinates of the alignment in both the reference and query sequences. For a given k-mer with start coordinates in both the reference and query sequences... , , The alignment, and the coordinates on the reference sequence are k-mer, query the k-mer coordinates that match it on the sequence. Calculated in the following way:

[0020] ,

[0021] in and The alignment lengths on the reference sequence and the query sequence are respectively, by and Calculations show that, in obtaining each pair of k-mers... and Then, the conflict score CS is calculated as follows:

[0022] .

[0023] In some embodiments, according to the method for filling gaps in genome assembly according to the present invention, the maximum alignment ratio PLA is calculated as the ratio of the longest alignment length to the ideal alignment region length.

[0024] In some embodiments, the method for filling gaps in genome assembly according to the present invention includes the calculation of the total length of conflict-free alignments (LNA), the proportion of conflict-free alignments (PNA), and the number of conflict-free alignment matching bases (MBN), comprising:

[0025] Overlapping alignments are merged, with the length of the overlapping region on both the reference and query sequences being greater than a specified length, and the difference between the overlap lengths on the reference and query sequences being less than a specified value. The merged alignment is then compared with the longest alignment to detect conflicting alignments.

[0026] For each alignment, calculate its distance d from the longest alignment on both the reference sequence and the query sequence. r and d q If max(d) r d q )>2min(d r d q If it is a conflict comparison, then it is considered a conflict comparison.

[0027] The distance between alignments is defined by subtracting the coordinate of the earlier-ending alignment from the coordinate of the later-starting alignment; after eliminating conflicting alignments, the total length of conflict-free alignments (LNA), the proportion of conflict-free alignments (PNA), and the number of matching bases in the conflict-free alignments (MBN) are calculated using the remaining conflict-free alignments.

[0028] In some embodiments, the method for filling gaps in genome assembly according to the present invention further includes the step of calculating a uniform reliability score URS, which is calculated based on two or more combinations of conflict score CS, maximum alignment ratio PLA, total length of conflict-free alignment LNA, proportion of conflict-free alignment PNA, and number of conflict-free matching bases MBN.

[0029] Preferably, the Unified Reliability Score (URS) is calculated using the following formula:

[0030] URS=β0+β1CS+β2PLA+β3LNA+β4PNA+β5MBN,

[0031] Where (β) is the weight parameter, which can be obtained through training with known data.

[0032] In some embodiments, the method for filling gaps in genome assembly according to the present invention further includes providing a reliability score kMAPQ based on a unique K-mer for each individual alignment, which is calculated by the following formula:

[0033]

[0034] Among them, l LIS The length of the longest increasing subsequence is represented by lr, the length of the reference sequence is represented by lq, and the length of the query sequence is represented by lq.

[0035] In some embodiments, the method for filling gaps in genome assembly according to the present invention, wherein the step of generating a set of candidate read sequences for filling gaps includes: constructing an overlap graph and performing a heuristic breadth-first search on the overlap graph to obtain a series of paths from one side of the gap region to the other side, each path representing a set of candidate read sequences for filling gaps;

[0036] Preferably, in the heuristic breadth-first search process, when a vertex v is visited, the algorithm randomly selects a fixed number of k unvisited subsequent vertices based on the probability calculated from the weights of the edges connected to it. The algorithm terminates when the search reaches the vertex representing the other side.

[0037] A second aspect of the invention provides a system for filling gaps in genome assembly, comprising a user end and a server end, wherein the following steps are performed on the server end:

[0038] (1) Provide chromosome-level assembly data and whole-genome third-generation sequencing data containing multiple read sequences generated from third-generation sequencing data;

[0039] (2) Generate a first alignment file and a second alignment file, wherein the first alignment file contains alignment data between reads and chromosomes, and the second alignment file contains alignment data between different reads;

[0040] (3) Filter out the read segments that do not have gap regions and retain the read segments with gap regions, wherein the gap regions include the gap region and the two sides of the gap;

[0041] (4) Process the first comparison file and the second comparison file to filter out conflicting comparison data and obtain valid comparison data;

[0042] (5) Evaluate the quality of each alignment in the valid alignment data obtained in step (4) and generate a reliability score;

[0043] (6) Generate a set of candidate read segments for filling gaps using a heuristic breadth-first search based on reliability scores (optional);

[0044] Specifically, the personal terminal is designed to read files from the server based on an index file construction mechanism. Because the alignment files are too large and slow to read, an index file is generated for each alignment file. Since the alignment files are ordered, for the first alignment file, the index file records the file position of the first occurrence of a chromosome, allowing direct location of the required chromosome during reading. For the second alignment file, the index file records the file positions of all sequencing sequence names, allowing direct location of the required read segment. The contents of the index file are then compiled into a binary file and transmitted to the PC for independent distribution and execution.

[0045] The method and system of this invention overcome the technical problems existing in traditional manual gap-filling processes. It allows for efficient gap-filling during genome assembly without the need for specialized technicians, requiring only simple operation, and is particularly suitable for telomere-to-telomere level genome assembly. Simultaneously, this invention significantly reduces gap-filling time, avoiding high workload. Correctness, time efficiency, and space efficiency were tested using an Arabidopsis thaliana reference genome. Attached Figure Description

[0046] Figure 1 A schematic diagram of an exemplary gap-filling process of the present invention.

[0047] Figure 2 Exemplary diagram illustrating conflict alignment. For a set of alignment results for a reference sequence and a query sequence, if the alignment start position of alignment A in the reference sequence is larger than the alignment start position of alignment B, while in the query sequence, the alignment start position of alignment A is smaller than the alignment start position of alignment B, this indicates a conflict between the two alignment results A and B. A and B are called a conflicting alignment.

[0048] Figure 3 The diagram illustrates the maximum alignment percentage. For a set of alignment results {a1, a2, ..., an} between the reference sequence and the query sequence, the alignment length of ai is the largest among the alignment results. ai is called the maximum alignment in this set of alignments. However, ai may not be an ideal alignment. The maximum alignment is used to obtain the ideal alignment position. For the attached diagram, the ideal alignment position on the reference sequence is from 1 to r. e +(q L -q e The position of q in the query sequence is q. s -r s to q L The maximum comparison ratio is shown in the following formula:

[0049]

[0050] Figure 4 For a set of alignment results for the reference sequence and the query sequence, regardless of whether it is the reference sequence or the query sequence, the alignment termination position of alignment A is larger than the alignment start position of alignment B, and the alignment termination position of alignment B is larger than the alignment termination position of alignment A. This indicates that there is an overlap between the two alignment results A and B. The two alignments A and B can be merged into one alignment. The starting position of the merged alignment is the starting position of A, and the ending position is the ending position of B. Detailed Implementation

[0051] Various exemplary embodiments of the present invention will now be described in detail. This detailed description should not be considered as a limitation of the present invention, but rather as a more detailed description of certain aspects, features, and embodiments of the present invention.

[0052] It should be understood that the terminology used in this invention is merely for describing particular embodiments and is not intended to limit the invention. Furthermore, with respect to numerical ranges in this invention, it should be understood that the upper and lower limits of the range and each intermediate value between them are specifically disclosed. Any stated value or intermediate value within a stated range, as well as each smaller range between any other stated value or intermediate value within said range, are also included in this invention. The upper and lower limits of these smaller ranges may be independently included or excluded from the range.

[0053] Unless otherwise stated, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. While only preferred methods and materials have been described herein, any methods and materials similar or equivalent to those described herein may be used in the implementation or testing of this invention. All references to this specification are incorporated by way of citation to disclose and describe methods and / or materials associated with those references. In the event of any conflict with any incorporated reference, the content of this specification shall prevail.

[0054] Example

[0055] Methods for filling gaps in genome assembly

[0056] In a first aspect, this invention develops a method to assist in the manual gap-filling process, saving time and effort during genome assembly. Even those without bioinformatics skills can efficiently perform manual gap-filling using, for example, a personal computer with just a mouse click, and directly obtain a complete assembly containing the filled sequences. The method of this invention only requires input of chromosome-level assemblies generated from an automated workflow, alignments between long-read sequencing sequences and the assemblies, and pairwise sequence alignments. For each gap to be filled, this invention provides two modes for gap-filling: (1) a "stepwise" mode and (2) an "automatic" mode.

[0057] "Gradual" model

[0058] In the "stepwise" mode, after obtaining a set of candidate read sequences, the method of the present invention allows a technician or general system operator to select a specific (or reliable) candidate read sequence from the set according to their needs, thereby obtaining a long-read sequencing sequence that can be extended. This fills one or both sides of the gap, resulting in an updated gap and both sides. The technician or general system operator performs the filling step continuously or repeatedly in the same manner until ligation is achieved on both sides.

[0059] In some embodiments, the present invention provides a set of standards and a uniform score for each read to assess the reliability of read alignments to one or both sides, preferably further generating a set of recommended read sequences for selection by technicians or general system operators. Such embodiments are advantageous for those with less knowledge of genome assembly and sequence alignment, particularly general system operators. In some preferred embodiments, although the recommended reliability is high, technicians or general system operators can still select any read for filling as needed.

[0060] In some embodiments, the system of the present invention is designed to have visualization capabilities, enabling any step or operation of the invention to be demonstrated graphically or animatedly. For example, the gap region and its surrounding region (or sides), as well as long read sequences comparing one or both sides, or reversal actions, can be visualized, allowing skilled personnel or general operators to visually understand and be aware of the assembly process or progress.

[0061] In some implementations, after a technician or general system operator selects one or more specific (or reliable) candidate read sequences from the candidate read sequence set according to their own needs and fills one or both sides of the gap, if it is found that the sequence filling cannot continue due to the selection of a false positive read in a previous step of the filling process, a step-by-step "undo" operation can be performed to return to the first step of the incorrect selection and reselect.

[0062] In some implementations, to facilitate filling larger gaps, the system of the present invention allows for interruption during the filling process and saving of the filling state, with the operation resuming the next time by reloading the state.

[0063] "Automatic" mode

[0064] In the "automatic" mode of this invention, the invention preferably automatically predicts at least one possible gap-filling solution. Each solution comprises a series of overlapping long read sequences. In the case of multiple solutions, the automatic mode will produce multiple sets of overlapping long read sequences. The most reliable list of read sequences can then be selected for filling, either automatically by machine or manually by an operator as needed. In automatic mode, after filling a gap, the sequence of the filled region can be automatically generated by the backend program and added to the chromosome-level assembly of the FASTA file.

[0065] Systems for filling gaps in genome assembly

[0066] Because automated genome assembly processes cannot generate complete chromosome sequences, almost all assemblies involve time-consuming manual gap-filling processes to fill in missing sequences from complex genomic regions (e.g., tandem and scattered repeats). In most assemblies, automated processes take only a few days, while manual processes can take months or even years. Furthermore, manual work requires advanced information skills. Addressing the significant conflict between the limited computing and storage capabilities of personal computers and the massive amounts of data (such as sequences and alignments) at the whole genome level, the system of this invention improves time and space efficiency through optimization. The system of this invention is designed to include two separate parts. The first part, the server-side (exemplarily based on a first operating system such as Linux), handles all time-consuming algorithmic steps for preprocessing, while the second part, the personal end (i.e., the PC, which can be based on a second operating system such as Windows), is designed to perform only simple file operations and other functions such as visualization. However, achieving this goal is not easy due to the dependencies between gap-filling steps. For example, during candidate sequence extension, the extension is based on alignments of reads with the new flanking sequences obtained in the previous step, which is very time-consuming. Since the newly obtained unknown sequences on both sides cannot be pre-calculated on the server side, the system of this invention employs a strategy that replaces the comparison between the read segments and the newly extended sides by comparing the read segments pairwise, thereby converting the repeated sequence comparisons on the PC side into a one-time comparison on the server side.

[0067] To improve space efficiency, the system of this invention employs a strategy on the server side to filter out reads and corresponding alignments in non-gap regions. Secondly, to evaluate the quality and reliability of candidate reads and their alignments with both sides, and to recommend reliable reads to system users or operators, this invention develops a k-mer-based method to assess the degree of conflict between segment alignments, based on which multiple scores can be calculated and displayed. Furthermore, by learning a linear regression model on a reference database, a weighted average of these multiple scores or combinations thereof can be calculated to obtain a unified score.

[0068] The system of this invention can select either the "automatic" mode or the "stepwise" mode described above. In the "automatic" mode, for example, to automatically generate candidate sequences of read segments for filling gaps, this invention constructs an overlap graph, where vertices represent read segments and their sides, edges exist between overlapping read segments, and a breadth-first search (BFS) algorithm is used to search for paths between the vertices representing the sides, where each path represents a candidate sequence, i.e., a combination of read segments. To address the problem that BFS cannot find all paths in a finite time for large gaps with numerous vertices and edges, this invention modifies the traditional BFS algorithm by adding a heuristic strategy. Edges corresponding to more reliable overlaps are visited with a higher probability, thus more reliable paths can be found earlier when executing the DFS algorithm.

[0069] Since the gaps during assembly are mostly repetitive genomic regions, accurate visualization of repetitive sequence alignment is crucial. Due to the high similarity of repetitive units, repetitive region alignment is often segmented, and the same sequence can be aligned with multiple regions. To remove erroneous alignments, this invention uses a specific algorithm as described above to find the optimal subset of non-conflicting alignments and remove other alignments.

[0070] In summary, a second aspect of the present invention provides a system for filling gaps in genome assembly, comprising a server side and a user side (sometimes abbreviated herein as "PC side"), wherein,

[0071] The following steps are performed on the server side:

[0072] (1) Provide chromosome-level assembly data and whole-genome third-generation sequencing data containing multiple read sequences generated from third-generation sequencing data;

[0073] (2) Generate a first alignment file and a second alignment file, wherein the first alignment file contains alignment data between reads and chromosomes, and the second alignment file contains alignment data between different reads;

[0074] (3) Filter out the read segments that do not have gap regions and retain the read segments with gap regions, wherein the gap regions include the gap region and the two sides of the gap;

[0075] (4) Process the first comparison file and the second comparison file to filter out conflicting comparison data and obtain valid comparison data;

[0076] (5) Evaluate the quality of each alignment in the valid alignment data obtained in step (4) and generate a reliability score;

[0077] (6) Generate a set of candidate reads for filling gaps using a heuristic breadth-first search based on reliability scores (optional).

[0078] Specifically, the personal terminal is designed to read files from the server based on an index file construction mechanism, compile the contents of the input file into an executable file, and transmit it to the server for independent distribution and execution.

[0079] In an exemplary embodiment, the PC client of the present invention is designed to have layout organization, style settings, system prompts, and user interaction functions. These functions can be implemented using methods known in the art, for example, by writing code in Python or by utilizing known toolkits, such as using packages like PyQt5, QFluentWidgets, and PyQtGraph to design a graphical user interface (GUI).

[0080] In some embodiments, the PC client of this invention is designed with iterative file traversal capabilities, enabling automatic reading of required files. To improve software runtime efficiency, the PC client of this invention is designed with an index file construction mechanism. Based on the content of the input file, images conforming to specific requirements, such as SVG format, are generated, and these images are subsequently automatically displayed using, for example, the QtSvg module in PyQt5. Exemplarily, by using, for example, the Nuitka package, the files are compiled into executable files suitable for, for example, a first operating system such as Windows, thereby enabling independent distribution and execution on a server. In an exemplary embodiment, by implementing remote graphical user interface (GUI) display technology, a GUI can be provided on a second operating system such as Linux, ensuring its normal operation and good usability in the second operating system environment.

[0081] In an exemplary embodiment, the personal terminal of this invention is designed for visualization and is configured with two basic modes: the first mode is mainly used to display the alignment between reads and chromosomes, which is the initial stage of constructing the assembly pathway; the second mode is used to display the alignment between reads, which is an extension of the assembly pathway. Based on the relative position of the alignment to the gap region, reads are further divided into three different display methods: at the left end of the gap, at the right end of the gap, and able to cross the gap. The operator can choose any one of these methods for subsequent selection. Simultaneously, reads are arranged in descending order of alignment reliability. When the operator clicks on a read, the alignment information related to that read is automatically displayed, and various parameters measuring alignment reliability are listed to intuitively assist the operator in selection. When determining candidate reads, the design can include two methods: one is to directly add the read to the assembly path; the other is to add it to a candidate list, so that different alignment situations are displayed simultaneously to facilitate selection.

[0082] Test case

[0083] Using genome assembly results and high-fidelity long sequencing reads provided in the National Genome Database PRJCA007112, tests were conducted on the model organism Arabidopsis thaliana. A 100 kbp gap was artificially created in the centromere region of Arabidopsis chromosome 1. All results below are based on this dataset.

[0084] 1. Filtering

[0085]

[0086] Regarding the filtering of reads from non-gap regions to improve time and space efficiency, the overlap file obtained before filtering was 795G, while after filtering, the overlap file was only 86G, reducing space usage by 90%. Furthermore, the accuracy of the filtering results was tested. To determine if reads truly belong to gap regions, reads were aligned to the Arabidopsis T2T genome, and alignments belonging to gap regions were extracted. Then, using the same conditions as filtering non-gap region reads—alignment quality greater than 10 and alignment length greater than 500bp—the alignments were filtered, ultimately yielding 1117 reads belonging to gap regions. In terms of read quantity, this test filtered 90% of non-gap region reads; in terms of accuracy, only 5 of the 1117 correct reads were filtered, retaining over 99% of the gap region reads. The final results indicate that filtering is very effective in reducing space usage and can retain the vast majority of gap region reads without significantly impacting subsequent gap filling.

[0087] 2. Recommended score

[0088] CS PLA LNA PNA MBM Correlation 9.74％ 12.46％ 12.64％ 19.02％ 34.98％

[0089] The recommendation score test was based on the alignment of reads with the complete Arabidopsis genome and the overlap between reads. Minimap2 was used for both read-to-genome alignment and overlap acquisition. For positive samples, this test utilized the alignment results between reads and the complete Arabidopsis genome. Alignments with a quality score greater than 30 and an alignment length greater than 500 bp were considered high-quality alignments. Furthermore, if two high-quality alignments corresponded to reads with an overlap length greater than 1000 bp in the complete genome, the overlap between these two reads was considered a positive sample and recorded as 1. For negative samples, reads that showed alignment information between two reads in the overlapping read files but were not included in the positive samples were considered problematic and recorded as 0. After selecting positive and negative samples, the correlation between the five provided recommendation scores and the results for both positive and negative samples was calculated. The final results show that even the conflict score with the lowest correlation has a 9.74% correlation with the correctness of overlap between read segments, while the highest number of matching bases has a 34.98% correlation. This indicates that the recommended scores provided by this invention can help users select appropriate read segments for extension when filling gaps.

[0090] Although the invention has been described with reference to exemplary embodiments, it should be understood that the invention is not limited to the disclosed exemplary embodiments. Various adjustments or changes may be made to the exemplary embodiments described in this specification without departing from the scope or spirit of the invention. The scope of the claims should be interpreted in the broadest possible sense to cover all modifications and equivalent structures and functions.

Claims

1. A method for filling gaps in genome assembly, characterized in that, Includes the following steps: (1) Provide chromosome-level assembly data and whole-genome third-generation sequencing data containing multiple read sequences generated from third-generation sequencing data; (2) Generate a first alignment file and a second alignment file, wherein the first alignment file contains alignment data between reads and chromosomes, and the second alignment file contains alignment data between different reads; (3) Filter out the read segments that do not have gap regions and retain the read segments with gap regions, wherein the gap regions include the gap region and the two sides of the gap; (4) Process the first alignment file and the second alignment file to filter and remove conflicting alignment data to obtain valid alignment data. The step of filtering and removing conflicting alignments includes: in the reference sequence and query sequence In the fragment comparison set A between the two, the comparisons will be made ( According to its position in the reference sequence The corresponding fragment on ( The reference coordinates are sorted in ascending order to obtain the query sequence. The corresponding segment ( ) sequence They are arranged in the same order, and simultaneously obtain reference coordinates and Same order ( Reference coordinate sequence of the sequence ;Calculate using binary search algorithm Longest increasing subsequence By collecting all from And it is in the query sequence The corresponding segment ( The reference coordinates are contained in The comparison within the subset yields the optimal conflict-free subset. This filters out conflict comparisons. The reference coordinates include start and / or end coordinates, and their longest increasing subsequence... The following condition must be met: for any two values ​​in the sequence and ,like exist Previously there was ; (5) Evaluate the quality of each alignment in the valid alignment data obtained in step (4) and generate a reliability score; (6) Input the chromosome-level assembly results, sequencing data, first and second alignment files and reliability score files into the PC, and select the appropriate sequence to fill the gaps.

2. The method for filling gaps in genome assembly according to claim 1, characterized in that, Further includes at least one of the following steps: In step (1), data for generating assembly units for chromosome-level assembly is provided; Based on reliability scores, a heuristic breadth-first search is used to generate a set of candidate read segments to fill the gaps.

3. The method for filling gaps in genome assembly according to claim 1, characterized in that... The quality of the alignment is assessed by calculating at least one of the following: Conflict Score (CS), Maximum Alignment Ratio (PLA), Total Collision-Free Alignment Length (LNA), Collision-Free Alignment Ratio (PNA), and Number of Matched Bases in Collision-Free Alignment (MBN). The conflict score CS is calculated as follows: A certain number of k-mers are randomly selected within the alignment region of the reference sequence, and their positions in the query sequence are calculated. The position of the k-mer in the query sequence is obtained based on the start and end coordinates of the alignment in both the reference and query sequences. For a given k-mer with start coordinates in both the reference and query sequences... , , The alignment, and the coordinates on the reference sequence are k-mer, query the k-mer coordinates that match it on the sequence. Calculated in the following way: ， in and The alignment lengths on the reference sequence and the query sequence are respectively, by and Calculations show that, in obtaining each pair of k-mers... and Then, the conflict score CS is calculated as follows: CS ; The maximum alignment ratio (PLA) is calculated as the ratio of the longest alignment length to the ideal alignment region length. The calculation of the total length of collision-free alignment (LNA), the proportion of collision-free alignment (PNA), and the number of collision-free matching bases (MBN) includes: Overlapping alignments are merged, with the length of the overlapping region on both the reference and query sequences being greater than a specified length, and the difference between the overlap lengths on the reference and query sequences being less than a specified value. The merged alignment is then compared with the longest alignment to detect conflicting alignments. For each alignment, calculate its distance from the longest alignment on both the reference sequence and the query sequence. and ,like If so, it is considered a conflict comparison; The distance between alignments is defined by subtracting the coordinate of the earlier-ending alignment from the coordinate of the later-starting alignment; after eliminating conflicting alignments, the total length of conflict-free alignments (LNA), the proportion of conflict-free alignments (PNA), and the number of matching bases in the conflict-free alignments (MBN) are calculated using the remaining conflict-free alignments. The method further includes the step of calculating the Uniform Reliability Score (URS), which is calculated based on two or more combinations of the following: Collision Score (CS), Maximum Alignment Ratio (PLA), Total Collision-Free Alignment Length (LNA), Collision-Free Alignment Ratio (PNA), and Number of Matched Bases in Collision-Free Alignment (MBN). The Unified Reliability Score (URS) is calculated using the following formula: ， in, to These are the weight parameters, which can be obtained through training with known data.

4. The method for filling gaps in genome assembly according to claim 1, characterized in that, This further includes providing a reliability score kMAPQ based on a unique K-mer for each individual comparison, which is calculated by the following formula: ， Among them, l LIS The length of the longest increasing subsequence is represented by lr, the length of the reference sequence is represented by lq, and the length of the query sequence is represented by lq.

5. The method for filling gaps in genome assembly according to claim 1, characterized in that, The steps for generating a set of candidate read sequences for filling the gap include: constructing an overlap graph and performing a heuristic breadth-first search on the overlap graph to obtain a series of paths from one side of the gap region to the other side, each path representing a set of candidate read sequences for filling the gap.

6. The method for filling gaps in genome assembly according to claim 5, characterized in that, In the heuristic breadth-first search process, when a vertex is visited... At that time, the algorithm calculates the probability based on the weights of the edges connected to it, and randomly selects a fixed number of... The algorithm terminates when the search reaches a vertex representing the other side, after finding an unvisited subsequent vertex.

7. A system for filling gaps in genome assembly, characterized in that, It includes a personal terminal and a server terminal, wherein the method for filling gaps in genome assembly as described in any one of claims 1-6 is performed on the server terminal.

Images