A reference genome based assembly method for single molecule real time sequencing

Abstract

The application relates to the field of individual genome assembly, in particular to a reference genome-based assembly method for single-molecule real-time sequencing. The application aims to solve the problems of high raw error rate of single-molecule real-time sequencing data, large calculation resource and memory consumption of traditional de novo assembly, difficulty in crossing repeat sequences and complex structure variation regions in a genome and the like. The process comprises the following steps: step one, constructing a hash index based on a reference genome; adopting a phased alignment strategy to construct an alignment skeleton of long read sequences; recording breakpoints of the alignment skeleton of the long read sequences as SV signals; step two, clustering the SV signals; step three, constructing an overlapping graph of the long read sequences; converting the overlapping graph into a directed string graph and simplifying the directed string graph; and replacing the reference sequence with long read sequence alignment information.

Description

Technical Field

[0001] This invention relates to the field of individual genome assembly, and more specifically to a reference genome-based assembly method for real-time single-molecule sequencing. Background Technology

[0002] The reference genome-based assembly strategy utilizes a high-quality reference genome already existing in the same individual as a guide, and sequence alignment assists in the assembly of the target genome. This strategy can effectively address challenges such as repetitive sequence regions and complex structural variations by leveraging known reference genome information, thereby significantly improving the efficiency and accuracy of individual genome assembly, while also enabling rapid identification of structural variations.

[0003] The rise of single-molecule real-time sequencing (represented by PacBio SMRT and ONT) has brought revolutionary breakthroughs to genomics research. Leveraging its long read length, this technology effectively overcomes the limitations of traditional short-read sequencing in resolving repetitive sequences and complex genomic variation regions. However, the early single-molecule real-time sequencing data had a raw error rate as high as approximately 15%, posing new challenges to genome assembly. De novo assembly based on single-molecule real-time sequencing requires extremely high computational resources and memory, prompting researchers to develop reference genome-guided assembly strategies. These strategies can anchor sequencing fragments to the correct positions in the reference genome, fully utilizing the known structure and location of repetitive or complex variation regions on the reference genome to guide assembly algorithms to correctly traverse repetitive regions and identify variation sites between different individuals, greatly improving the quality, continuity, and integrity of genome assembly.

[0004] Single-molecule real-time sequencing, with its ability to generate data instantly, is suitable for reference genome-guided assembly, which is fast and computationally inefficient. It can dynamically assess the quality and integrity of the assembly results during sequencing, and terminate the sequencing process once the target of genome assembly and variant detection is detected. This efficient mode of simultaneous sequencing and assembly significantly reduces sequencing costs, providing a cost-effective solution for individual genome assembly. Summary of the Invention

[0005] The purpose of this invention is to address the problems of high raw error rate in single-molecule real-time sequencing data, large computational and memory consumption of traditional de novo assembly, and difficulty in traversing repetitive sequences and complex structural variation regions in the genome. Therefore, this invention proposes a reference genome-based assembly method for single-molecule real-time sequencing.

[0006] The specific process of a reference genome-based assembly method for single-molecule real-time sequencing is as follows:

[0007] Step 1: Construct a hash index based on the reference genome; construct the alignment backbone of long read sequences using a phased alignment strategy; record the alignment backbone breakpoints of long read sequences as SV signals; the specific process is as follows:

[0008] Step 11: Construct a hash index based on the reference genome;

[0009] Steps 1 and 2: Based on the hash index, construct the alignment skeleton of the long read sequence using a phased alignment strategy:

[0010] Phase 1: Construct the initial candidate alignment skeleton and output preliminary alignment results in real time;

[0011] Phase 2: Optimize the alignment skeleton of long read sequences by re-recalling duplicate minimizers;

[0012] Step 13: Record the alignment skeleton breakpoints of the long read sequence as SV signals;

[0013] Step 2: Cluster the SV signals; the specific process is as follows:

[0014] Step 2: Construct a candidate window for SV signals;

[0015] Step 22: Based on the candidate window of the SV signal, merge adjacent SV windows;

[0016] Step 3: Construct a long read sequence overlap graph; convert the overlap graph into a directed string graph and simplify it; replace the reference sequence with the long read sequence alignment information; the specific process is as follows:

[0017] Step 3: 1. Construct a long read sequence overlap graph;

[0018] Step 32: Convert the overlapping graph into a directed string graph and simplify it;

[0019] Step 33: Based on the simplified directed string graph from Step 32, replace the reference sequence with the long read sequence alignment information.

[0020] The beneficial effects of this invention are as follows:

[0021] This invention constructs a genome index and aligns sequences. This step proposes a two-stage, efficient sequence alignment method. By constructing a hash index based on the reference genome and progressively optimizing the alignment backbone, it significantly reduces computational complexity while maintaining accuracy. In a real-time data stream environment, an efficient hash index is first constructed based on the reference genome. Then, a staged alignment strategy is adopted: the first stage focuses on constructing an initial candidate alignment backbone and outputs preliminary alignment results in real time for subsequent analysis; the second stage optimizes the alignment backbone by re-recalling duplicate minimizers (k-mers with the smallest hash value selected through a sliding window), achieving progressively refined mapping. This strategy, through a progressively refined alignment process, effectively avoids the computational overhead of traditional full-sequence alignment, enabling faster and more accurate sequence alignment results even in real-time scenarios.

[0022] This invention utilizes long-read sequence clustering based on SV breakpoints. This step employs local clustering to identify candidate variant windows, avoiding global genome alignment analysis and thus reducing computational resource consumption. First, the reference genome is divided into fixed-size windows, and based on real-time acquired long-read sequence alignment information, detected SV signals are immediately assigned to the corresponding windows. Subsequently, candidate windows containing significant SV signals are selected, and an adaptive window fusion strategy is used to dynamically integrate adjacent candidate windows based on the spatial distribution characteristics of the variant signals. This strategy, through local focusing and dynamic screening, effectively avoids the computational burden of traditional whole-genome scanning, enabling more efficient and accurate identification of structural variants.

[0023] This invention is based on local assembly-based sequence reconstruction. This step employs a local assembly strategy to generate high-quality individual genome sequences by progressively detecting and integrating overlapping information. During real-time assembly, candidate sequences within a window are first compared to each other to detect overlap information in real time. Then, near structural variation signals, local sequence assembly is performed based on the detected overlap information, resulting in continuous and accurate local sequences, with real-time feedback on assembly progress. Finally, using long-read alignment information, the high-quality locally assembled sequences replace corresponding regions in the reference genome, achieving correction and integration of the target genome. This progressive mapping and real-time correction approach not only improves the continuity of assembly but also better preserves structural variation information between individuals.

[0024] This invention provides sequencing data with high accuracy, low computational cost, low memory consumption, and the ability to traverse repetitive sequences and complex structural variation regions in the genome. It solves the problems of high raw error rate in single-molecule real-time sequencing data, as well as the high computational and memory consumption of traditional de novo assembly and its difficulty in traversing repetitive sequences and complex structural variation regions in the genome. Attached Figure Description

[0025] Figure 1 This is a flowchart of the present invention;

[0026] Figure 2 This is a schematic diagram of the process of the present invention;

[0027] Figure 3 To construct a directed string graph;

[0028] Figure 4 This is a transformation graph of directed strings, where a, b, c, d, e, and f represent long-read sequence numbers. Detailed Implementation

[0029] Specific implementation method one: Combining Figure 1 , Figure 2 This embodiment describes a reference genome-based assembly method for single-molecule real-time sequencing. The specific process is as follows:

[0030] This invention proposes a highly efficient assembly method based on a reference genome for single-molecule real-time sequencing. By integrating long-read sequencing data with reference genome information, it achieves accurate reconstruction of an individual genome. In real-time sequencing scenarios, this method can progressively construct the assembly map during data generation and output partial data in real time for immediate analysis and feedback. It mainly includes the following three core steps:

[0031] Step 1: Construct a hash index based on the reference genome; construct the alignment backbone of long read sequences using a phased alignment strategy; record the alignment backbone breakpoints of long read sequences as SV signals;

[0032] The specific process is as follows:

[0033] This invention proposes an efficient method for aligning third-generation sequencing data with a reference genome. In real-time data processing scenarios, it achieves immediate data feedback and dynamic optimization through phased output of partial data and progressive mapping analysis. This method utilizes a minimizer to construct a hash index of the reference genome, enabling rapid location of alignment regions. A two-round alignment strategy is employed to build the alignment backbone, significantly improving alignment accuracy while substantially reducing computational resource consumption. In real-time processing mode, the system can output partial intermediate results in real time during data generation and dynamically adjust subsequent alignment processes based on online analysis. By finely integrating local matching blocks and utilizing the recall mechanism of the repetition minimizer, this method not only effectively addresses complex repetitive regions and structural variations but also possesses high flexibility and scalability.

[0034] Step 11: Construct a hash index based on the reference genome;

[0035] Steps 1 and 2: Based on the hash index, construct the alignment skeleton of the long read sequence using a phased alignment strategy:

[0036] Phase 1: Construct the initial candidate alignment skeleton and output preliminary alignment results in real time;

[0037] The second stage: optimize the alignment skeleton of long read sequences by re-recalling duplicate minimumizers (k-mers with the smallest hash value selected by a sliding window);

[0038] Step 13: Record the alignment skeleton breakpoints of the long read sequence as SV signals;

[0039] Step 2: Cluster the SV signals; the specific process is as follows:

[0040] A local clustering approach is employed to identify candidate variant windows, avoiding global genome alignment analysis and thus reducing computational resource consumption. First, the reference genome is divided into fixed-size windows, and based on real-time long-read sequence alignment information, detected SV signals are immediately assigned to the corresponding windows. Then, candidate windows containing significant SV signals are selected, and an adaptive window fusion strategy is used to dynamically integrate adjacent candidate windows based on the spatial distribution characteristics of the variant signals. This strategy, through local focusing and dynamic selection, effectively avoids the computational burden of traditional whole-genome scanning, enabling more efficient and accurate identification of structural variants.

[0041] Step 2: Construct a candidate window for SV signals;

[0042] Step 22: Based on the candidate window of the SV signal, merge adjacent SV windows;

[0043] Step 3: Construct a long read sequence overlap graph; convert the overlap graph into a directed string graph and simplify it; replace the reference sequence with the long read sequence alignment information; the specific process is as follows:

[0044] In real-time data processing scenarios, the system not only progressively constructs local overlap maps and assembly maps, but also outputs some intermediate results in real time during the assembly process. Its workflow includes three main steps: constructing a long-read sequence overlap map for a given window, performing local sequence assembly around SV signals, and replacing sequences in the reference genome with long-read sequence alignment information. This workflow fully leverages the advantages of long-read data spanning complex repetitive regions and structural variations, improving the continuity and accuracy of assembly while reducing computational complexity, and supporting real-time data analysis and feedback.

[0045] Step 3: 1. Construct a long read sequence overlap graph;

[0046] Step 32: Convert the overlapping graph into a directed string graph and simplify it;

[0047] Step 33: Based on the simplified directed string graph from Step 32, replace the reference sequence with the long read sequence alignment information;

[0048] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that step one involves constructing a hash index based on a reference genome; the specific process is as follows:

[0049] A minimizer-based hash index was constructed for the reference genome;

[0050] 1) Given a reference genome, a window size of w (default value is 15bp), and a step size of 1 base;

[0051] The reference genome is divided into windows starting from the first base of the reference genome, resulting in all windows;

[0052] 2) Calculate w k-mers within each window using a hash function. The k-mer with the smallest hash value is defined as the minimumizer for the corresponding window. Collect all minimumizers from the reference genome.

[0053] 3) Sort all minimizers of the reference genome from front to back according to base position and add them to the linear list L, and construct a hash table H;

[0054] The key of hash table H is l-mer (default value: 11bp), and the value of hash table H is a list of positions of the minimizer in the linear list L.

[0055] Through the above steps, this invention constructs an efficient two-level hash index structure, providing a foundation for fast searching in subsequent comparisons.

[0056] The other steps and parameters are the same as in Specific Implementation Method 1.

[0057] Specific Implementation Method Three: This implementation method differs from Specific Implementation Method One or Two in that, in steps one and two, a phased comparison strategy is used based on hash indexes to construct the comparison skeleton of long read sequences:

[0058] Phase 1: Construct the initial candidate alignment skeleton and output preliminary alignment results in real time;

[0059] The second stage: optimize the alignment skeleton of long read sequences by re-recalling duplicate minimumizers (k-mers with the smallest hash value selected by a sliding window);

[0060] The specific process is as follows:

[0061] Sequence alignment based on the alignment backbone employs a two-round alignment strategy to efficiently align the input long read sequence with the reference genome, while supporting real-time data output and stepwise mapping analysis. First, an initial alignment backbone is constructed using minimizer hash matching and sparse dynamic programming. Then, a repeated minimizer recall mechanism is used to optimize the alignment results, significantly improving the alignment accuracy of complex genomic regions while maintaining computational efficiency.

[0062] 1) Construct the initial candidate alignment skeleton and output the preliminary alignment results in real time; the specific process is as follows:

[0063] 11) For a given long read sequence Q, a window size of w (default value is 15bp) and a step size of 1 base;

[0064] Divide the long read sequence Q into windows to obtain all windows;

[0065] When constructing the hash table of the reference genome, if a minimumr appears many times (>R) in the reference, it is likely located in a repetitive region, leading to numerous matching positions and potentially interfering with backbone construction. To avoid this interference, when extracting minimumrs from candidate sequences and matching them with the reference hash table, the number of matching positions in the corresponding reference is checked. If a minimumr appears more than R times in the reference genome, the window corresponding to that minimumr is disregarded.

[0066] 12) Calculate w k-mers in each window using a hash function. The k-mer with the smallest hash value is defined as the minimumizer of the corresponding window. Collect the minimumizers of all windows in the long read sequence Q.

[0067] 13) Match the minimizer of each window of the long read sequence Q with the hash index of the reference genome to locate the candidate alignment position of the minimizer of each window in the long read sequence Q in the reference genome;

[0068] Using the same window size w and hash function as when constructing the reference genome index, minimizers are extracted from Q and matched against the hash index of the reference genome; the candidate alignment positions of the minimizers in Q in the reference genome are efficiently located through fast lookup of the hash index.

[0069] 14) Let the i-th matching pair after matching in 13) be (q i ,r i );

[0070] Where, q iThis represents the position of the long read sequence Q corresponding to the i-th matching pair; r i This indicates the position of the reference genome corresponding to the i-th matching pair;

[0071] If two adjacent matching pairs (q) i ,r i ) and (q i+n ,r i+n ) satisfies |(q i+n -q i )-(r i+n -r i If |≤ε′, then the pair (q) is considered a match. i ,r i ) and (q i+n ,r i+n ) have collinearity, matching pairs (q) i ,r i ) and (q i+n ,r i+n The matching blocks (MBs) are merged to form larger matching blocks, where MB = {(q,r)|q∈[q...}. start ,q end ], r∈[r start ,r end ]};

[0072] Where ε′ represents the allowable relative position deviation threshold, that is, the maximum error that can be tolerated in the difference between the displacements of two matching pairs in the query sequence and the reference sequence;

[0073] q i+n r represents the position of the long read sequence Q corresponding to the (i+1)th matching pair. i+n This indicates the position of the reference genome corresponding to the (i+1)th matching pair;

[0074] i represents the index of the current matching pair in all extracted minimizer matching pairs in 3), and n represents the matching pair interval used for collinearity testing;

[0075] q represents the position of the matching block in Q, and r represents the position of the matching block in the reference genome; q start Indicates the starting position of q, q end Indicates the ending position of q, r start Indicates the starting position of r, r end Indicates the ending position of r;

[0076] These matching blocks represent regions where long read sequences are highly consistent with the reference genome.

[0077] 15) Construct a directed acyclic graph (DAG), treating each matching block (MB) as a node in the DAG. If two matching blocks (MB) are contiguous on the reference genome, and the order in which they appear on the long read sequence Q is consistent with their order on the reference genome (q on the long read sequence Q), then the matching blocks are considered a node in the DAG. i In q i+n Previously, on the reference genome r i In r i+n (Previously), then in the i-th matching block MB i up to the j-th matching block MB j Establish a directed edge between them and define the edge weight w. ij and penalty weight p ij ;

[0078] The following recursive formula is calculated sequentially for each matching block according to the coordinate order on the reference genome, MB. i Before, MB j Following this; the initial score for the first matching block is 0;

[0079] The score for each MB is calculated using sparse dynamic programming with the following recursive formula:

[0080]

[0081] in,

[0082] S(MB j ) represents the j-th matching block MB j The score;

[0083] S(MB i ) represents the i-th matching block MB i The score;

[0084] Predecessors represents all connections that can be directly linked to MB. i Front-end matching block (in MB) i The preceding set of matching blocks;

[0085] Typically, a candidate sequence will match multiple regions in the reference genome.

[0086] 16) Select the matching block with the highest score as the last point of the path, according to the above formula. Proceed backwards until the first matching block is found, obtain the path, and use the path as a candidate alignment skeleton;

[0087] The candidate alignment skeleton is transformed into an alignment between the candidate sequence and the reference sequence. The alignment information includes the candidate sequence identifier, the position of the candidate sequence in the reference sequence, and the alignment skeleton score (the score of the matching block corresponding to the last point on the path).

[0088] The candidate sequence is the currently aligned long read sequence;

[0089] 2) Optimize the candidate alignment skeleton; the specific process is as follows:

[0090] After constructing candidate alignment skeletons in the first stage, the main goal of the second stage is to optimize the alignment skeleton by re-recalling duplicate minimizers and to detect possible structural variations. In real-time scenarios, the system will monitor and dynamically analyze key breakpoints in the preliminary alignment skeleton in real time.

[0091] 21) In order to reduce the interference of repetitive region comparison on the results, if a minimumer appears more than R times in the reference genome, the minimumer that appears more than R times will not be considered.

[0092] 22) Because long read sequences in repetitive regions may lack effective matching information, this can lead to false positive breakpoints in the alignment backbone. To address this issue, in 16) the candidate alignment backbones output in real time retain matching blocks whose alignment backbone scores exceed a preset threshold;

[0093] 23) For any gap between two adjacent matching blocks, recall all minimizers in the gap (those duplicate minimizers that were previously filtered out because they appeared more than R times);

[0094] In the first stage, to avoid noise caused by high repetition, minimizers that appear more than R times in the reference genome were actively ignored. However, this would lead to a lack of matching information in some regions (especially between anchor points), resulting in false positive breakpoints in the alignment backbone. Therefore, in the second stage, these previously ignored minimizers need to be reactivated.

[0095] 24) Merge the corrected candidate alignment skeletons.

[0096] Other steps and parameters are the same as in specific implementation method one or two.

[0097] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that the merged and corrected candidate alignment skeleton in step 24) is as follows:

[0098] 241) Determine whether the re-recalled minimizer and the matching block preceding the gap interval are collinear in the reference genome and length sequence;

[0099] 242) Determine whether the re-recalled minimizer and the matching block following the gap interval are collinear in the reference genome and length sequence;

[0100] 243) If 2231) and 2232) are collinear, the minimumizer that is recalled is used as the new matching block; execute 2234);

[0101] Otherwise, the minimumizer that is recalled again will not be used as a new matching block;

[0102] 244) Merge the matching blocks retained in the candidate alignment skeleton in 22) with the new matching blocks to generate an ordered list of matching blocks;

[0103] 245) Determine whether the ordered matching block list generated in 244) needs to be merged to obtain the final matching block list, which serves as the alignment skeleton after merging and correction; the specific process is as follows:

[0104] 2451) Read adjacent matching block pairs in the ordered matching block list generated in 244);

[0105] 2452) When the positions of adjacent matching blocks in the reference genome and in the long read sequence overlap, or the gap length between adjacent matching blocks is less than the preset merging threshold Δ, the adjacent matching blocks will be merged into one matching block.

[0106] The starting point of the merged matching block in the reference genome is the minimum starting coordinate of the adjacent matching block in the reference genome; the starting point of the merged long read sequence is the minimum starting coordinate of the adjacent matching block in the long read sequence; the ending point of the merged matching block in the reference genome is the maximum ending coordinate of the adjacent matching block in the reference genome; and the ending point of the merged long read sequence is the maximum ending coordinate of the adjacent matching block in the long read sequence.

[0107] 2453), until all adjacent matching blocks have been traversed, output the merged and corrected alignment skeleton.

[0108] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0109] Specific Implementation Method 5: This implementation method differs from one of the specific implementation methods 1 to 4 in that the skeleton breakpoints are recorded as SV signals in step 13.

[0110] The specific process is as follows:

[0111] Determine if there are any breakpoints in the optimized skeleton;

[0112] If there is a breakpoint, the breakpoint is used as a candidate SV signal, and the position information of the candidate SV signal is recorded (R). s ,R e G s G eSimultaneously record the split alignment information (start and end positions of the breakpoint on the reference genome, start and end positions of the breakpoint on the long read sequence, forward and reverse strands of the long read sequence, start and end positions of the long read sequence, and start and end positions of the long read sequence on the reference genome).

[0113] Among them, R s and R e G represents the start and end positions of the long read sequence, respectively. s and G e These represent the start and end positions on the reference genome, respectively;

[0114] If there are no breakpoints, record the full-length alignment information (start and end positions on the reference genome, start and end positions of your own genome, long read sequence forward and reverse strands, alignment backbone score);

[0115] Full-length alignment information refers to a complete skeleton obtained after the above-mentioned supplementary matching points. It has no breaks, the skeleton can be connected from beginning to end in a consistent direction, and it covers the entire long read. This complete skeleton information is the full-length alignment information.

[0116] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0117] Specific Implementation Method Six: This implementation method differs from Specific Implementation Methods One to Five in that the candidate window for the SV signal is constructed in step Two-One;

[0118] The specific process is as follows:

[0119] 1) Divide the reference genome into a series of overlapping windows of length Lbp (default: 1024bp), where the overlap length of each window is L / 2;

[0120] In real-time data processing scenarios, as long-read sequencing data is continuously input, the system dynamically distributes the detected SV signals (long-read sequencing data compared with the reference genome, which is backbone breakpoint information) to various windows;

[0121] 2) Select candidate windows containing sufficient SV signals by using either condition one or condition two;

[0122] Condition 1: For a certain window Θ of 1), if the following conditions are met... Then window Θ is considered a candidate window;

[0123] Among them, C sv RD represents the number of SV signals within window Θ. Θ P represents the long read sequence coverage depth of window Θ. Θ This is a custom parameter (default value: 6), typically representing the sequencing depth. Ensure that, at greater coverage depths, the number of SV signals significantly exceeds the background noise;

[0124] Condition 2: For a certain window Θ in 1), due to the repetitive nature of the genome, sequencing errors, or coverage gaps, RD Θ The value of can be very low or even zero. In this case, window Θ is considered a candidate window as long as it contains an SV signal. This condition ensures that potential SV signals are not missed in low-coverage areas.

[0125] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0126] Specific Implementation Method Seven: This implementation method differs from one of Specific Implementation Methods One to Six in that, in step two, the candidate windows based on the SV signal are merged into adjacent SV windows;

[0127] The specific process is as follows:

[0128] To further optimize the reliability of candidate windows, this invention merges windows surrounding multiple adjacent SV signals. A series of candidate windows exist on the reference genome, consisting of SV signals and their surrounding regions. Each candidate window for an SV signal is denoted as Q. k =[R s,k ,R e,k ], k = 1, 2, ..., K;

[0129] Among them, R s,k and R e,k These represent the start and end positions of the k-th candidate window on the reference genome, respectively;

[0130] Each candidate window Q k For one or more (two or more) long read sequences, each candidate window Q k The identifiers (each long read sequence within each window has one identifier) ​​constitute a set M. k ;

[0131] To optimize window reliability, this invention employs a greedy strategy to merge adjacent windows. The merging conditions include:

[0132] Adjacency of positions: The spacing between adjacent windows satisfies Δ k =R s,k+1 -R e,k When ≤δ, the adjacent window Q k With Q k+1 Merge into a new window Q′ k =[R s,k ,R e,k+1 ];

[0133] Where δ is the distance threshold, indicating that the two windows are sufficiently close; Δ k This represents the gap between the starting position of the (k+1)th window and the ending position of the kth window;

[0134] R s,k+1 R represents the starting coordinates of the (k+1)th window on the reference genome. e,k This represents the end coordinates of the k-th window on the reference genome;

[0135] Long read sequence connectivity: Adjacent windows are connected by the same long read sequence, i.e. Connect adjacent windows Q k With Q k+1 Merge into a new window Q′ k =[R s,k ,R e,k+1 ], and update each window Q k The set of identifiers is M′=M k ∪M k+1 .

[0136] By integrating and dynamically merging neighboring windows supported by the same long read sequence in real time, a series of high-confidence local windows are obtained to integrate the scattered SV signals. At the same time, the merging results are output in real time for further real-time analysis and dynamic adjustment.

[0137] The other steps and parameters are the same as those in one of the specific implementation methods one to six.

[0138] Specific Implementation Method Eight: This implementation method differs from Specific Implementation Methods One through Seven in that step three-one involves constructing a long-read sequence overlap graph; the specific process is as follows:

[0139] Based on the set of long read sequences in the remaining window of step 22, an overlapping graph is constructed using the string graph method;

[0140] In the initial string diagram (such as) Figure 4 In this process, each non-containment sequencing sequence is split into two vertices—fB and fE (f represents the sequencing sequence, fB represents the start position of the sequencing sequence, and fE represents the end position of the sequencing sequence), corresponding to the 5′ and 3′ positions of the read, respectively. There is no need to explicitly store the overall orientation or complementary sequence of the read; instead, the Watson-Crick complement for the read is determined during subsequent assembly based on the vertex type (B or E). Specifically, if a path passes through vertex v (v.type = B), the reverse complement of v.read is taken; if v.type = E, the forward sequence is preserved (see...). Figure 4 The dual-vertex model can accurately represent the location of both ends of a segment and flexibly handle the splicing relationship of positive or negative chains. Figure 4The genome sequence shown therein contains two key regions: the regions marked in yellow represent unique sequences in the genome, each of which occurs only once in the entire genome, and the regions marked in red and green represent repeat sequences, wherein sequences of the same color occur twice in the genome and belong to identical repeat units.

[0141] For any overlap o between two reads f and g, if the corresponding overlap interval on f is [o.f.beg,o.f.end] and o.f.beg<o.f.end, two directed edges need to be added in the string graph to represent the overhang sequences outside the overlap:

[0142] 1. Left extension edge: it starts from vertex f.B and points to g.B, with an edge label of f[0,o.f.beg], which is the remaining base sequence before the overlap region of f;

[0143] 2. Right extension edge: it starts from vertex g.E and points to f.E, with an edge label of comp(g[o.g.end,g.len]), wherein comp(...) represents Watson-Crick complementation, that is, the remaining base sequence after the overlap region of g is complemented.

[0144] When the overlap is in the opposite direction on g (o.g.beg>o.g.end), the same rule is also applied, except that when selecting the non-overlap interval, it is necessary to first locate the corresponding segment according to the overlap direction. In this way, whether the overlap is forward or reverse, each edge accurately marks how two reads are spliced outside the overlap, ensuring that when walking along any path, the overlap region can be automatically aligned, and the remaining overhang sequences and the reads themselves can be spliced in order according to the correct direction and strand information;

[0145] The present invention adopts an improved strategy to construct a long read overlap graph, which is based on minimizer technology and can identify overlap regions between long reads. In order to avoid missing short SV signals, strict restriction conditions are imposed when constructing an alignment skeleton;

[0146] The overlap graph satisfies the following constraints:

[0147] Sequence overlap length constraint: for overlapping fragments of two long reads within one window, the coordinates on the long reads satisfy q k+1,start -q k,end ≤Δ max ;

[0148] wherein, Δ m ax is a preset maximum allowable interval, q k+1,start represents the starting position of the overlapping fragment of the (k+1)-th long read; q k,end represents the end position of the overlapping fragment of the k-th long read;

[0149] This constraint ensures that adjacent matching blocks are close enough in length sequence to merge consecutive matching information and obtain more accurate information.

[0150] Reference genome position overlap constraint: Two long read sequences within a single window must have overlapping coordinate ranges on the reference genome, and the distance between the coordinates of the two overlapping long read sequences on the reference genome must satisfy |r... k+1,start -r k,end |≤D max ;

[0151] Among them, D max This is the preset maximum distance to avoid incorrect connections caused by excessively long distances between matching blocks, ensuring consistency between the aligned skeleton and the real reference structure; k+1,start This represents the starting coordinates of the (k+1)th long read overlap segment on the reference genome; r k,end This indicates the end coordinates of the overlapping segment of the k-th long read on the reference genome;

[0152] In real-time scenarios, this module can gradually output local overlap information and perform statistical analysis to obtain overlap data close to the base level.

[0153] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0154] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One to Eight in that, in step three two, the overlapping graph is converted into a directed string graph and simplified.

[0155] The specific process is as follows:

[0156] In the second step, this invention converts the long-read sequence overlap graph into a directed string graph and generates contigs derived from linear single paths in the graph through a series of graph simplification operations. The real-time constructed graph allows the system to gradually output partial graph structures during the assembly process. The operations in this stage mainly include several key steps: First, by removing transitive edges, redundant information in the graph is eliminated, simplifying the structural complexity of the graph; second, short branches caused by sequencing errors are pruned to simplify the graph structure.

[0157] First, transitive edges are removed to eliminate redundancy and simplify the complexity of the graph. This invention employs a transitive reduction algorithm to simplify directed string graphs, aiming to simplify the graph representation of long-read sequences by precisely removing redundant connections.

[0158] 1) Perform edge deletion on the overlapping graph obtained in step 31; the specific process is as follows:

[0159] All nodes (the length sequence includes 1 positive chain and 1 negative chain, corresponding to 2 nodes) are initialized to an unvisited state to ensure that each node is visited only once during processing;

[0160] Let E be the set of outgoing edges of a node u in an overlapping graph (the set of all directed edges originating from node u and pointing to other nodes). out (u)={e1,e2,…,e q ,…,e n};

[0161] Among them, each edge e q The length is l(e) q );

[0162] Calculate the length of the longest outgoing edge And set a fault tolerance parameter For each edge e q If there exists a path P from u to v (u and v are edges e) q The two vertices satisfy Then it is believed that e q It is a redundant transitive edge, so edge e q By marking and removing redundant connections, the complexity of strings can be simplified. This process can gradually output the simplified graph structure during real-time graph building.

[0163] Where length(P) represents the length of path P;

[0164] 2) Prune short branches (tips) caused by sequencing errors in the overlapping graph obtained in 1) to improve the graph's accuracy. This algorithm is used to prune short tips in the genome assembly graph, which are usually caused by sequencing errors. The algorithm removes these potentially misleading structures by identifying and deleting short paths connected to a single ingress or egress edge. Short branches caused by sequencing errors indicate that the sequencing itself is flawed and may contain errors, resulting in short branches (short paths originating from an endpoint but unable to continue) in the graph. The specific process is as follows:

[0165] Let the short branch tip caused by sequencing errors be a path from endpoint v1→v2→…v q →…→v n The edge sequence {e1,e2,…,e q ,…,e n The total length of the tip is

[0166] If the tip length satisfies L tip <L min If the tip is not removed, it will reduce error messages and optimize the subsequent map interpretation process. The real-time output map data can help monitor the tip pruning effect.

[0167] If the tip length does not meet L tip <L min If so, the tip will not be deleted;

[0168] Among them, L min This represents the preset threshold; v1 represents the starting node of the tip, v q This represents the first intermediate node in the tip path, v n This indicates the termination node of the tip;

[0169] 3) Obtain the bubble structure based on the overlapping graph obtained in 2), and then reduce the size of the overlapping graph obtained in 2) based on the bubble structure;

[0170] The specific process is as follows:

[0171] 31) The assembly graph is a directed graph G = (V, E), where V is the set of nodes and E is the set of edges;

[0172] 32) Use depth-first search to traverse graph G. For each node v∈V, if the out-degree of node v is greater than 1, then v is considered a potential bubble entrance; otherwise, v is not considered a potential bubble entrance.

[0173] 33) For two candidate nodes u, v ∈ V, if there exist at least two disjoint paths P1 and P2 from u to v, and satisfy d(u, v) = min{length(P1), length(P2)} <D max If the structure of path P1 and P2 is true, then the structure of path P1 and P2 is defined as bubble structure B(u,v); otherwise, the structure of path P1 and P2 is not bubble structure B(u,v).

[0174] Where d(u,v) represents the shortest path length from u to v, D max The preset maximum detection distance; length(P1) represents the length of path P1; length(P2) represents the length of path P2;

[0175] 34) After detecting the bubble structure B(u,v), the algorithm enters the bubble popping stage, further distinguishing between two types based on the internal path of the bubble:

[0176] The main bubble path is denoted as P. main =argmaxS(P);

[0177] Where S(P) represents the number of SV signals contained in all paths P in the bubble structure B(u,v), P main Represents the major haplotype of the genome;

[0178] The alternative bubble path is denoted as P. alter S(P) alter ) <S(Pmain ), P alter Represents a minor haplotype;

[0179] Among them, S(P main ) indicates the main bubble path P main The number of SV signals included;

[0180] 4) Obtain a continuous linear contig.

[0181] The specific process is as follows:

[0182] 41) Generate a substitution sequence based on the substitution bubble path, delete the substitution sequence in the overlapping graph obtained in 3), and obtain a new overlapping graph;

[0183] This invention ensures the accuracy of the splicing path and extracts the globally optimal splicing path through the following methods:

[0184] 42) Check whether all nodes on the remaining paths in the new overlap graph obtained in 41) satisfy the condition that they all originate from the same continuous segment of the same reference genome, and that the relative directions (forward and backward strands) of two adjacent nodes appearing on the long read sequence are consistent with the relative directions (forward and backward strands) of their appearance on the reference genome.

[0185] If satisfied, a continuous linear contig is obtained;

[0186] If not satisfied, then execute 43);

[0187] 43) In the new overlapping graph obtained in 41), identify the source node (connecting the source node to other nodes) with an in-degree of zero and an out-degree of greater than zero, and the sink node (connecting other nodes to the sink node) with an in-degree (number of incoming edges) of greater than zero and an out-degree (number of outgoing edges) of zero. Use Dijkstra's algorithm to calculate the shortest path from the source node to the sink node, thus obtaining a continuous linear contig; achieving the globally optimal contig in extreme repetitive and complex branching environments.

[0188] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0189] Specific Implementation Method 10: This implementation method differs from one of the specific implementation methods one to nine in that, in step three-three, the long read sequence alignment information is replaced with the reference sequence based on the simplified directed string graph in step three-two.

[0190] The specific process is as follows:

[0191] Replace the reference sequence with the contig based on the start and end positions of the contig in the reference genome;

[0192] Determine the termination conditions:

[0193] (1) Convergence of mutation detection: Repeat steps two to three. When the number of new mutations in the data is less than 5 for three consecutive times, it indicates that the main mutation area has been fully covered.

[0194] (2) Replacement coverage stability: When the genome sequence coverage of the replaced region reaches 95% and the replacement fluctuation range (the number of replaced reference sequences as a percentage of the total number of reference sequences) is less than 2% when steps two to three are repeated twice, it indicates that the key region update has become stable.

[0195] If any condition is met, the real-time sequencing process will automatically terminate, and the final updated sequence along with all substitution sequences will be integrated and stored in a new FASTA file for output. This dual dynamic termination mechanism ensures the efficiency of the analysis process and the reliability of the results.

[0196] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0197] Example:

[0198] Experimental setup

[0199] This invention utilizes both real and simulated datasets. The real dataset used is the PacBio HiFi dataset, derived from the widely studied GIAB sample HG002 (NA24385). The data was then downsampled using seqtk (5x, 10x, 20x, 30x), for a total of five real datasets. The simulated dataset used SURVIVOR to randomly simulate 675 and 732 structural variants, targeting both Arabidopsis thaliana and Drosophila melanogaster. The variant length ranged from 50 to 10000 bp, the homozygosity ratio of the structural variants was set to 0.6, and the coverage was 5x, 10x, 20x, 30x, and 50x.

[0200] Table 1. Genomic information of species and information of generated datasets

[0201]

[0202] This invention outputs individual genome sequences (FASTA format), then uses minimap2 to align the output sequences to a reference genome (parameter "-ax asm10–cs-r2k"). Next, samtools generates sorted and indexed BAM files, and SVIM-asm is used to call SV in haploid mode to output VCF files. Finally, bcftools is used for sorting, bgzip compression, and tabix indexing. The SV detection performance of all five tools is evaluated using Truvari (parameters "-r 1000-p 0--dup-to-ins--pickmulti--passonly").

[0203] When evaluating the predictive performance of the model, we use metrics including precision, recall, and F1 score. These metrics comprehensively evaluate the detection effect by quantifying the model's ability to distinguish between true variants (positive examples) and non-variant sites (negative examples).

[0204] 1. Accuracy

[0205] Accuracy measures the proportion of actual variants among the variant sites predicted by the model. This metric directly reflects the reliability of the variant detection results.

[0206]

[0207] TP represents the true variant sites correctly predicted by the model, while FP represents the normal sites incorrectly predicted by the model. High accuracy means that the variant sites reported by the model have high reliability.

[0208] 2. Recall rate

[0209] Recall represents the proportion of samples that are actually positive that are correctly predicted as positive by the model. Recall assesses the model's ability to detect true variants.

[0210]

[0211] FN represents the true variant sites that the model failed to detect.

[0212] 3. F1 value

[0213] The F1 score is the harmonic mean of precision and recall, which comprehensively reflects the model's performance in terms of both accuracy and coverage. When there is a trade-off between precision and recall, it can provide a unified performance metric.

[0214]

[0215] When both accuracy and recall are high, the F1 score will also be high; however, if one of them is low, the F1 score will be significantly affected, reflecting an imbalance in model performance.

[0216] 4. Experimental Results

[0217] This invention selected four state-of-the-art structural variation (SV) detection tools (cuteSV2 v2.1.1, sniffles2 v2.2, SVDSS v2.0.0, and SVIM v2.0.0) for comparison, and the experimental results are summarized in the table below.

[0218] Table 2 shows the results: accuracy, recall, and F1 score.

[0219]

[0220]

[0221] This invention may have other embodiments. Without departing from the spirit and essence of this invention, those skilled in the art can make various corresponding changes and modifications according to this invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A reference genome-based assembly method for single-molecule real-time sequencing, characterized in that: The specific process of the method is as follows: Step 1: Construct a hash index based on the reference genome; construct the alignment backbone of long read sequences using a phased alignment strategy; record the alignment backbone breakpoints of long read sequences as SV signals; the specific process is as follows: Step 11: Construct a hash index based on the reference genome; Steps 1 and 2: Based on the hash index, construct the alignment skeleton of the long read sequence using a phased alignment strategy: Phase 1: Construct the initial candidate alignment skeleton and output preliminary alignment results in real time; Phase 2: Optimize the alignment skeleton of long read sequences by re-recalling duplicate minimizers; Step 13: Record the alignment skeleton breakpoints of the long read sequence as SV signals; Step 2: Cluster the SV signals; the specific process is as follows: Step 2: Construct a candidate window for SV signals; Step 22: Based on the candidate window of the SV signal, merge adjacent SV windows; Step 3: Construct a long read sequence overlap graph; convert the overlap graph into a directed string graph and simplify it; replace the reference sequence with the long read sequence alignment information; the specific process is as follows: Step 3:

1. Construct a long read sequence overlap graph; Step 32: Convert the overlapping graph into a directed string graph and simplify it; Step 33: Based on the simplified directed string graph from Step 32, replace the reference sequence with the long read sequence alignment information; The specific process of steps one and two is as follows: Step 121: Construct the initial candidate alignment skeleton and output the preliminary alignment results in real time; the specific process is as follows: Step 1211: For a given long read sequence Given a window size of w and a step size of 1 base; For long read sequences Divide the window into all windows; Step 1212: Calculate w k-mers within each window using a hash function. The k-mer with the smallest hash value is defined as the minimumizer for the corresponding window. Collect long read sequences. Minimizer for all windows; Steps 1-3: Transform the long read sequence The minimizer of each window is matched against the hash index of the reference genome to locate long read sequences. Candidate alignment positions of the minimizer in the reference genome for each window; Steps 1-4: Let the i-th matching pair after matching in steps 1-3 be... ; in, This represents the long read sequence corresponding to the i-th matching pair. Location; This indicates the position of the reference genome corresponding to the i-th matching pair; If two adjacent matching pairs and satisfy Then it is considered a match. and Having collinearity will match pairs and Merge to form a matching block, matching block ; in, This indicates the permissible relative position deviation threshold; This represents the long read sequence corresponding to the (i+1)th matching pair. Location, This indicates the position of the reference genome corresponding to the (i+1)th matching pair; This indicates the sequence number of the current matching pair among all the minimizer matching pairs extracted in steps one through three. This represents the matching pair interval used for collinearity testing; express The position of the matching block express The position of the matching block; express The starting position, express The ending position, express The starting position, express The ending position; Step 1215: Construct a directed acyclic graph and connect each matching block. Consider it as a node in a directed acyclic graph, if two matching blocks The positions on the reference genome are contiguous, and the two matching blocks are... If the order of occurrence of the long read sequence Q is consistent with the order of occurrence of its occurrence on the reference genome, then in the i-th matching block... up to the j-th matching block Establish a directed edge between them and define the edge weight. and penalty weight ; Calculate the score for each MB: in, This represents the j-th matching block. The score; Represents the i-th matching block The score; This indicates all that can be directly connected. The set of predecessor matching blocks; Step 1216: Select the highest-scoring matching block as the last point of the path. Derive backwards from the above formula until the first matching block is found to obtain the path. Use the path as the candidate alignment skeleton. The candidate alignment skeleton is transformed into an alignment between the candidate sequence and the reference sequence. The alignment information includes the candidate sequence identifier, the position of the candidate sequence in the reference sequence, and the alignment skeleton score information. The candidate sequence is the currently aligned long read sequence; Step 122: Optimize the candidate alignment skeleton; the specific process is as follows: Step 1: If a minimumr appears more than R times in the reference genome, then the minimumr that appears more than R times is not considered. Step 1222: In the candidate alignment skeleton output in real time in Step 1216, retain the matching blocks whose alignment skeleton scores exceed the preset threshold; Step 1223: For any gap between two adjacent matching blocks, recall all minimizers within the gap. Steps 1-4: Merge and correct the candidate alignment skeletons.

2. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 1, characterized in that: In each step, a hash index is constructed based on the reference genome; The specific process is as follows: Step 111: Given a reference genome, a window size of w, and a step size of 1 base; The reference genome is divided into windows starting from the first base of the reference genome, resulting in all windows; Step 112: Calculate w k-mers within each window using a hash function. The k-mer with the smallest hash value is defined as the minimumizer for the corresponding window. Collect all minimumizers of the reference genome. Step 113: Sort all minimizers of the reference genome from front to back according to base position and add them to the linear list L, and construct a hash table H; The key of hash table H is l-mer, and the value of hash table H is a list of positions of the minimizer in the linear list L.

3. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 2, characterized in that: In steps one through four, the candidate alignment skeletons are merged and corrected. The specific process is as follows: Step 12241: Determine whether the re-recalled minimizer and the matching block preceding the gap interval are collinear in the reference genome and long reads; Step 12242: Determine whether the re-recalled minimizer and the matching block following the gap interval are collinear on the reference genome and long read sequences; Step 12243: If Step 12241 and Step 12242 are collinear, the minimumizer that is recalled again will be used as the new matching block; then proceed to Step 12244. Otherwise, the minimumizer that is recalled again will not be used as a new matching block; Step 12244: Merge the matching blocks retained in the candidate alignment skeleton from Step 1222 with the new matching blocks to generate an ordered list of matching blocks; Step 1-5: Determine whether the ordered matching block list generated in Step 1-4 needs to be merged to obtain the final matching block list, which serves as the alignment skeleton after merging and correction. The specific process is as follows: Step 1, 2, 2, 4, 5: Read adjacent matching block pairs sequentially from the ordered matching block list generated in Step 1, 2, 2, 4, 5. Step 122452: When the positions of adjacent matching blocks in the reference genome and in the long read sequence overlap, or the gap length between adjacent matching blocks is less than the preset merging threshold Δ, the adjacent matching blocks are merged into 1 matching block. The starting point of the merged matching block in the reference genome is the minimum starting coordinate of the adjacent matching block in the reference genome; the starting point of the merged long read sequence is the minimum starting coordinate of the adjacent matching block in the long read sequence; the ending point of the merged matching block in the reference genome is the maximum ending coordinate of the adjacent matching block in the reference genome; and the ending point of the merged long read sequence is the maximum ending coordinate of the adjacent matching block in the long read sequence. Steps 1-3: Continue until all adjacent matching blocks have been traversed, and output the merged and corrected alignment skeleton.

4. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 3, characterized in that: In steps one and three, the skeleton breakpoints are recorded as SV signals; The specific process is as follows: Determine if there are any breakpoints in the optimized skeleton; If there is a breakpoint, the breakpoint is used as a candidate SV signal, and the position information of the candidate SV signal is recorded. Simultaneously record split alignment information; in, and These represent the start and end positions of the long read sequence, respectively. and These represent the start and end positions on the reference genome, respectively; If there are no breakpoints, record the full-length comparison information.

5. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 4, characterized in that: In step two, a candidate window for the SV signal is constructed. The specific process is as follows: Step 211: Divide the reference genome into overlapping windows of length Lbp, where the overlap length of each window is L / 2; As long-read sequencing data continues to be input, the detected SV signals are dynamically allocated to each window; Step 2: Select candidate windows containing SV signals using either condition 1 or condition 2. Condition 1: For a certain window in step 2.11 If satisfied Then the window Considered a candidate window; in, Display window The number of SV signals inside, Display window Long read sequence coverage depth, These are custom parameters; Condition 2: For a certain window in step 2.11 As long as the window It contains SV signals, which will open the window. It is considered a candidate window.

6. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 5, characterized in that: In step two, candidate windows based on SV signals are merged into adjacent SV windows. The specific process is as follows: Each candidate window in the candidate window of the SV signal is represented as ; in, and They represent the first The start and end positions of each candidate window on the reference genome; Each candidate window For one or more long read sequences, each candidate window The identifiers constitute a set ; A greedy strategy is used to merge adjacent windows. The conditions for merging include: Adjacency of positions: The spacing between adjacent windows satisfies At that time, adjacent windows and Merge into a new window ; in, Distance threshold; Indicates the first The starting position of the first window and the second The gap between the end positions of each window; Indicates the first The starting coordinates of each window on the reference genome. Indicates the first The end coordinates of each window on the reference genome; Long read sequence connectivity: Adjacent windows are connected by the same long read sequence, i.e. , adjacent windows and Merge into a new window and update each window The set of identifiers is .

7. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 6, characterized in that: In step three, an overlap graph of long read sequences is constructed. The specific process is as follows: Based on the set of long read sequences in the remaining window of step 22, an overlapping graph is constructed using the string graph method; Overlapping graphs satisfy the following constraints: Sequence overlap length constraint: This requires that the overlapping segment of two long read sequences within a window must lie on the coordinates of the long read sequence. ; in, It is the preset maximum allowable interval; This indicates the starting position of the overlapping segment of the (k+1)th long read sequence; This indicates the end position of the overlapping segment of the k-th long read sequence; Reference genome position overlap constraint: Two long read sequences within a window must have overlapping coordinate ranges on the reference genome, and the distance between the coordinates of the two overlapping long read sequences on the reference genome must satisfy... ; in, It is the preset maximum distance; This indicates the starting coordinates of the (k+1)th long read overlap segment on the reference genome; This represents the end coordinates of the overlapping fragment of the k-th long read sequence on the reference genome.

8. The method for assembling a reference genome for single-molecule real-time sequencing according to claim 7, characterized in that: In step three, the overlapping graph is converted into a directed string graph and then simplified. The specific process is as follows: 1) Perform edge deletion on the overlapping graph obtained in step 31; The specific process is as follows: Initialize all nodes to an unvisited state; Let the set of outgoing edges of a node u in an overlapping graph be . ; Among them, each edge The length is ; Calculate the length of the longest outgoing edge And set a fault tolerance parameter. For each edge If there exists a path P from u to v satisfying Then it is believed It is a redundant transitive edge, which will... Mark and remove; in, Indicates the length of path P; 2) Prune the short branch tips caused by sequencing errors in the overlap map obtained in 1). The specific process is as follows: Let the short branch tip caused by sequencing errors be a branch from the endpoint edge sequence The total length of the tip is ; If the tip length meets the requirements < Then delete the tip; If the tip length is not met < If so, the tip will not be deleted; in, Indicates a preset threshold; This indicates the starting node of the tip. This indicates the first intermediate node in the tip path. This indicates the termination node of the tip; 3) Obtain the bubble structure based on the overlapping graph obtained in 2), and then reduce the size of the overlapping graph obtained in 2) based on the bubble structure; The specific process is as follows: 31) The assembly diagram is a directed graph. , where V is the set of nodes and E is the set of edges; 32) Use depth-first search to traverse graph G, for each node If the out-degree of node v is greater than 1, then v is considered a potential bubble inlet; otherwise, v is not considered a potential bubble inlet. 33) For two candidate nodes If there exist at least two disjoint paths from u to v and And satisfy Then define the path and The structure is a bubble structure. Otherwise, the path and The structure is not a bubble structure. ; in, This represents the shortest path length from u to v; The preset maximum detection range; Representing a path length; Representing a path length; 34) Upon detecting bubble structure Then, based on the internal path of the bubble, two categories were further distinguished: The main bubble path is denoted as ; in, Represents bubble structure All paths Number of SV signals included; The alternative bubble path is denoted as ; in, Indicates the main bubble path The number of SV signals included; 4) Obtain a continuous linear contig; The specific process is as follows: 41) Generate a substitution sequence based on the substitution bubble path, delete the substitution sequence in the overlapping graph obtained in 3), and obtain a new overlapping graph; 42) Check whether all nodes on the remaining paths in the new overlapping graph obtained in 41) satisfy the condition that they all originate from the same continuous segment of the same reference genome, and that the relative direction of two adjacent nodes on the long read sequence is consistent with the relative direction of their appearance on the reference genome. If satisfied, a continuous linear contig is obtained; If not satisfied, then execute 43); 43) In the new overlapping graph obtained in 41), identify the source node with zero in-degree and greater than zero out-degree and the sink node with greater than zero in-degree and zero out-degree. Use Dijkstra's algorithm to calculate the shortest path from the source node to the sink node, thus obtaining a continuous linear contig.

9. The method for assembly based on a reference genome for single-molecule real-time sequencing according to claim 8, characterized in that: In step 33, the long read sequence alignment information is replaced with the reference sequence based on the simplified directed string graph from step 32. The specific process is as follows: Replace the reference sequence with the contig based on the start and end positions of the contig in the reference genome; Determine the termination conditions: (1) Convergence of mutation detection; (2) Stable replacement coverage; If any condition is met, the real-time sequencing process will be automatically terminated, and the final updated sequence and all replacement sequences will be integrated and stored in a new FASTA file for output.

Images