Single molecule real-time sequencing-oriented assembling method based on reference genome

Abstract

The invention discloses a single molecule real-time sequencing-oriented assembling method based on a reference genome, and relates to the field of individual genome assembling, in particular to a single molecule real-time sequencing-oriented assembling method based on the reference genome. The method aims at solving the problems that single-molecule real-time sequencing data is high in original error rate, traditional de novo assembly is large in computing resource and memory consumption, and repeated sequences and complex structure variation regions in genomes are difficult to cross. The method comprises the following steps: step 1, constructing a Hash index based on a reference genome; constructing an alignment skeleton of the long read sequence by adopting a staged alignment strategy; recording the comparison skeleton breakpoint of the long read sequence as an SV signal; step 2, clustering is carried out on the SV signals; step 3, constructing a long read sequence overlapping graph; converting the overlapped graph into a directed character string graph and simplifying the directed character string graph; and replacing the reference sequence with the long read sequence comparison information.

Description

Technical Field

[0001] The present invention relates to the field of individual genome assembly, and in particular to a reference genome-based assembly method for single-molecule real-time sequencing. Background Art

[0002] Reference genome-based assembly strategies utilize existing high-quality reference genomes from individuals of the same species as a guide, aiding in the assembly of the target genome through sequence alignment. This strategy leverages known reference genome information to effectively address challenges such as repetitive sequence regions and complex structural variations, significantly improving individual genome assembly efficiency and sequence accuracy while 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 genomic research. With its long sequence read length, this technology effectively overcomes the limitations of traditional short-read sequencing in resolving repetitive sequences and complex genomic variation regions. However, the raw error rate of early single-molecule real-time sequencing data was as high as approximately 15%, which also brought new challenges to genome assembly. De novo assembly based on single-molecule real-time sequencing requires extremely high computing resources and memory, prompting researchers to develop reference genome-guided assembly strategies. This strategy can anchor the sequencing fragments to the correct position of the reference genome, fully utilizing the known structure and position of the repetitive sequences or complex variation regions on the reference genome, guiding the assembly algorithm to correctly cross the repetitive regions and identify variation sites between different individuals, greatly improving the quality, continuity and completeness of the genome assembly.

[0004] Single-molecule real-time sequencing, with its instant data generation, is suitable for fast, computationally efficient reference genome-guided assembly. It dynamically assesses the quality and completeness of the assembly during sequencing, terminating the sequencing process once the genome assembly and variant detection targets are complete. This efficient, simultaneous sequencing and assembly process significantly reduces sequencing costs and provides a cost-effective solution for individual genome assembly. Summary of the Invention

[0005] The purpose of the present invention is to solve the problems of high original error rate of single-molecule real-time sequencing data, large computing resources and memory consumption of traditional de novo assembly, and difficulty in crossing repetitive sequences and complex structural variation regions in the genome, and to propose a reference genome-based assembly method for single-molecule real-time sequencing.

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

[0007] Step 1: Build a hash index based on the reference genome; use a phased alignment strategy to build a long-read sequence alignment skeleton; record the long-read sequence alignment skeleton breakpoints as SV signals; the specific process is as follows:

[0008] Step 1: Build a hash index based on the reference genome;

[0009] Step 1 and 2: Based on the hash index, a phased alignment strategy is used to build the alignment skeleton of the long read sequence:

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

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

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

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

[0014] Step 21: Construct candidate windows for SV signals;

[0015] Step 22: Based on the candidate windows 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 31: 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 in step 32, replace the reference sequence with the long read sequence alignment information.

[0020] The beneficial effects of the present invention are:

[0021] The present invention constructs a genome index and an alignment sequence; this step proposes a two-stage efficient sequence alignment method, which significantly reduces the computational complexity while ensuring accuracy by constructing a hash index of a reference genome and gradually optimizing the alignment skeleton strategy. In a real-time data stream environment, an efficient hash index is first constructed based on the reference genome. Then a phased alignment strategy is adopted: the first phase takes the construction of an initial candidate alignment skeleton as the main goal, and outputs the preliminary alignment results in real time for subsequent analysis; the second phase optimizes the alignment skeleton by re-recalling repeated minimizers (k-mers with the minimum hash value selected by a sliding window) to achieve gradually refined mapping. This strategy effectively avoids the computational overhead of traditional full sequence alignment through a gradually refined alignment process, so that sequence alignment results can be obtained more quickly and accurately in real-time scenarios.

[0022] The present invention is based on the long-read sequence clustering of SV breakpoints. This step uses a local clustering method to identify candidate variation windows, avoiding the global alignment analysis of the genome, thereby reducing the consumption of computing resources. First, the reference genome is divided into windows of fixed size, and based on the long-read sequence alignment information collected in real time, the detected SV signals are immediately assigned to the corresponding windows. Subsequently, candidate windows containing significant SV signals are screened out, and an adaptive window fusion strategy is adopted to dynamically integrate adjacent candidate windows based on the spatial distribution characteristics of the variation signals. This strategy effectively avoids the computational burden of traditional whole-genome scanning through local focusing and dynamic screening, and can identify structural variations more efficiently and accurately.

[0023] The present invention is based on sequence reconstruction based on local assembly. This step adopts a strategy based on local assembly to generate high-quality individual genome sequences by gradually detecting and integrating overlapping information. In the real-time splicing process, the candidate sequences in the window are first compared with each other, and the overlapping information between sequences is detected in real time. Subsequently, near the structural variation signal, local sequence assembly is performed based on the detected overlapping information to splice out continuous and accurate local sequences, and the assembly progress is fed back in real time. Finally, the high-quality sequence obtained by local assembly is replaced with the corresponding region in the reference genome with the help of long-read comparison information to achieve correction and integration of the target genome. Through this method of step-by-step mapping and real-time correction, not only the continuity of the assembly is improved, but also the structural variation information between individuals is better preserved.

[0024] The sequencing data of the present invention has high accuracy, low computational complexity, low memory consumption, and can span repetitive sequences and complex structural variation regions in the genome. It solves the problems of high original error rate of single-molecule real-time sequencing data, and large computational resources and memory consumption of traditional de novo assembly, and difficulty in spanning repetitive sequences and complex structural variation regions in the genome. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 Flowchart of the present invention;

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

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

[0028] Figure 4 is a conversion graph of directed strings, where a, b, c, d, e, and f represent long-read sequence numbers. DETAILED DESCRIPTION

[0029] Specific implementation method 1: Combination 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 paper proposes an efficient reference genome-based assembly method for single-molecule real-time sequencing. By integrating long-read sequencing data with reference genome information, it achieves accurate reconstruction of individual genomes. In the real-time sequencing scenario, this method can gradually build an assembly map during the data generation process and output partial data in real time for immediate analysis and feedback. It mainly includes the following three core steps:

[0031] Step 1: Build a hash index based on the reference genome; use a phased alignment strategy to build a long-read sequence alignment skeleton; record the alignment skeleton breakpoints of the long-read sequence as SV signals;

[0032] The specific process is:

[0033] The present invention proposes an efficient method for aligning third-generation sequencing data with a reference genome. In a real-time data processing scenario, by outputting partial data in stages and gradually building a map for analysis, instant feedback and dynamic optimization of the data are achieved. The method uses a minimizer to construct a hash index of the reference genome, which can quickly locate the alignment area, and adopts a two-round alignment strategy to construct an alignment skeleton, thereby significantly reducing computing resource consumption while significantly improving the accuracy of the alignment. In real-time processing mode, the system can output some intermediate results in real time during the data generation process, and dynamically adjust the subsequent alignment process in combination with online analysis. By meticulously integrating local matching blocks and utilizing the recall mechanism of the repeat minimizer, this method can not only effectively solve the problems of complex repetitive regions and structural variations, but also has high flexibility and scalability.

[0034] Step 1: Build a hash index based on the reference genome;

[0035] Step 1 and 2: Based on the hash index, a phased alignment strategy is used to build the alignment skeleton of the long read sequence:

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

[0037] Phase 2: Optimize the alignment backbone of the long read sequence by re-calling the repeated minimizer (k-mer with the minimum hash value selected by the sliding window);

[0038] Step 13: Record the alignment backbone 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 method is used to identify candidate variant windows, avoiding global alignment analysis of the genome and thus reducing the consumption of computing resources. First, the reference genome is divided into windows of fixed size, and based on the long-read sequence alignment information collected in real time, the detected SV signals are immediately assigned to the corresponding windows. Subsequently, candidate windows containing significant SV signals are screened out, 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 effectively avoids the computational burden of traditional whole-genome scanning through local focusing and dynamic screening, and can identify structural variations more efficiently and accurately.

[0041] Step 21: Construct candidate windows for SV signals;

[0042] Step 22: Based on the candidate windows 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 gradually constructs local overlap maps and assembly maps, but also outputs some intermediate results in real time during the assembly process. The process 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 the help of long-read sequence alignment information. This process fully utilizes the advantages of long-read data spanning complex repetitive regions and structural variations, improving the continuity and accuracy of the assembly while reducing computational complexity, while supporting real-time data analysis and feedback.

[0045] Step 31: 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 of step 32, replace the reference sequence with the long read sequence alignment information;

[0048] Specific embodiment 2: This embodiment differs from specific embodiment 1 in that a hash index is constructed based on a reference genome in step 1; 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 (the default value is 15bp) and a step size of 1 base;

[0051] Divide the reference genome into windows starting from the first base of the reference genome to obtain all windows;

[0052] 2) Use the hash function to calculate the w k-mers in each window. The k-mer with the smallest hash value is defined as the minimizer of the corresponding window. Collect all minimizers of 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 build 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 the position list of minimizer in linear list L.

[0055] Through the above steps, the present invention constructs an efficient secondary hash index structure, which provides a basis for rapid search for subsequent comparisons.

[0056] Other steps and parameters are the same as those in the first embodiment.

[0057] Specific embodiment three: This embodiment differs from specific embodiments one or two in that, in steps one and two, a phased comparison strategy is used based on a hash index to construct a comparison skeleton for long read sequences:

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

[0059] Phase 2: Optimize the alignment backbone of the long read sequence by re-calling the repeated minimizer (k-mer with the minimum hash value selected by the sliding window);

[0060] The specific process is:

[0061] The backbone-based sequence alignment uses a two-round alignment strategy to efficiently align input long-read sequences against a reference genome, while supporting real-time data output and step-by-step mapping analysis. First, an initial alignment backbone is constructed using minimizer hash matching and sparse dynamic programming. Subsequently, a repeated minimizer recall mechanism is employed to optimize the alignment results, significantly improving alignment accuracy in complex genomic regions while maintaining computational efficiency.

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

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

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

[0065] When constructing the hash table for the reference genome, if a minimizer appears a significant number of times (>R) in the reference, it is likely located in a repetitive region, which can lead to a high number of matching positions and potentially interfere with backbone construction. To mitigate this interference, when extracting a minimizer from a candidate sequence and matching it with the reference hash table, the number of matching positions in the corresponding reference is checked. If a minimizer appears >R times in the reference genome, the window corresponding to that minimizer is disregarded.

[0066] 12) Use the hash function to calculate the w k-mers in each window. The k-mer with the smallest hash value is defined as the minimizer of the corresponding window. Collect the minimizers of all windows of 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, and 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, extract minimizers from Q and match these minimizers with the hash index of the reference genome; through fast lookup of the hash index, efficiently locate the candidate alignment positions of the minimizers in Q in the reference genome;

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

[0070] Among them, q irepresents the position of the long read sequence Q corresponding to the i-th matching pair; r i 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 )|≤ε′, then the matching pair (q i ,r i ) and (q i+n ,r i+n ) has collinearity, matching pairs (q i ,r i ) and (q i+n ,r i+n ) are merged into larger matching blocks (MatchBlocks, MBs), matching blocks MB={(q,r)|q∈[q start ,q end ],r∈[r start ,r end ]};

[0072] Where ε′ represents the allowed relative position deviation threshold, that is, the maximum error that can be tolerated between the displacement differences of two matching pairs on the query sequence and the reference sequence;

[0073] q i+n represents the position of the long read sequence Q corresponding to the i+1th matching pair, r i+n Indicates the position of the reference genome corresponding to the i+1th matching pair;

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

[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 end position of q, r start Indicates the starting position of r, r end Indicates the end position of r;

[0076] These matching blocks represent regions with high consistency between the long-read sequences and the reference genome.

[0077] 15) Construct a directed acyclic graph and regard each matching block MB as a node in the directed acyclic graph. If the positions of two matching blocks MB on the reference genome are continuous and the order in which the two matching blocks MB appear on the long read sequence Q is consistent with the order in which they appear on the reference genome (q on the long read sequence Q). i In q i+n Previously, on the reference genome i In r i+n Before), then in the i-th matching block MB i To the jth matching block MB j Create a directed edge between them and define the edge weight w ij and penalty weight p ij ;

[0078] Calculate the following recursive formula for each matching block in the order of coordinates on the reference genome: MB i In front, MB j Afterwards; the initial score of the first matching block is 0;

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

[0080]

[0081] in,

[0082] S(MB j ) represents the jth matching block MB j score;

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

[0084] Predecessors represent all the nodes that can be directly connected to the MB i Predecessor matching block (row in MB i The preceding matching block) collection;

[0085] Often 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. Derived forward until the first matching block is found, the path is obtained, and the path is used as the candidate alignment skeleton;

[0087] The candidate alignment skeleton is converted 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, the alignment skeleton score (the score of the matching block corresponding to the last point on the path), and other information;

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

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

[0090] After the candidate alignment skeleton is constructed in the first phase, the main goal of the second phase is to optimize the alignment skeleton by recalling repeated minimizers and detecting 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 repeated region comparison results, if a minimizer appears more than R times in the reference genome, the minimizer that appears more than R times will not be considered;

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

[0093] 23) For the gap interval between any two adjacent matching blocks, recall all minimizers in the gap interval (those duplicate minimizers that were previously filtered out because the number of occurrences exceeded R);

[0094] In the first stage, to avoid noise caused by high repetitions, minimizers with a number of occurrences > R in the reference genome are actively ignored. However, this will result in a lack of matching information in some regions (especially between anchor points), forming false positive breakpoints in the alignment backbone. Therefore, in the second stage, these previously ignored minimizers should be re-enabled.

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

[0096] Other steps and parameters are the same as those in the first or second embodiment.

[0097] Specific embodiment 4: This embodiment differs from one of specific embodiments 1 to 3 in that the corrected candidate alignment skeletons are merged in 24); the specific process is:

[0098] 241), determine whether the recalled minimizer and the matching block in front of the gap interval maintain collinearity on the reference genome and length sequence;

[0099] 242), determine whether the recalled minimizer and the matching block after the gap interval maintain collinearity on the reference genome and length sequence;

[0100] 243), when 2231) and 2232) maintain collinearity at the same time, the recalled minimizer is used as a new matching block; execute 2234);

[0101] Otherwise, the recalled minimizer is not used as a new matching block;

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

[0103] 245) Determine whether the ordered matching block list generated in 244) needs to be merged to obtain the final matching block list as the merged and corrected comparison skeleton; the specific process is:

[0104] 2451), sequentially reading 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 the positions in the long read sequence overlap, or the length of the gap between adjacent matching blocks is less than the preset merging threshold Δ, the adjacent matching blocks are merged into one matching block;

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

[0107] 2453) until all adjacent matching blocks are traversed and the merged and corrected comparison skeleton is output.

[0108] The other steps and parameters are the same as those in the first to third embodiments.

[0109] Specific embodiment 5: This embodiment differs from any one of specific embodiments 1 to 4 in that in step 13, the skeleton breakpoints are recorded as SV signals;

[0110] The specific process is:

[0111] Determine whether there are breakpoints in the alignment skeleton after optimization;

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

[0113] Among them, R s and R e Represent the starting and ending positions on the long read sequence, G s and G e Represent the starting and ending positions on the reference genome, respectively;

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

[0115] Full-length alignment information refers to the complete skeleton obtained after the above-mentioned supplementary matching points, without any breakpoints. The skeleton can be connected from beginning to end in a consistent direction and covers the entire long read. The information of this complete skeleton is the full-length alignment information.

[0116] The other steps and parameters are the same as those in the first to fourth embodiments.

[0117] Specific embodiment 6: This embodiment differs from any one of specific embodiments 1 to 5 in that the candidate window of the SV signal is constructed in step 21;

[0118] The specific process is:

[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 the real-time data processing scenario, the system dynamically distributes the detected SV signals (the comparison of long-read sequencing data with the reference genome, that is, the backbone breakpoint information) to each window as long-read sequencing data continues to be input;

[0121] 2) Select candidate windows containing sufficient SV signals through condition 1 or condition 2;

[0122] Condition 1: For a window Θ in 1), if Then window Θ is regarded as a candidate window;

[0123] Among them, C sv Indicates the number of SV signals within the window Θ, RD Θ represents the long read sequence coverage depth of window Θ, P Θ It is a custom parameter (default value: 6), which is generally the depth of sequencing. Ensure that the number of SV signals can significantly exceed the background noise at a high coverage depth;

[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, as long as the window Θ contains SV signals, the window Θ is considered a candidate window. 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 the first to fifth embodiments.

[0126] Specific embodiment seven: This embodiment differs from any one of specific embodiments one to six in that, in step 22, adjacent SV windows are merged based on the candidate windows of the SV signal;

[0127] The specific process is:

[0128] In order to further optimize the reliability of the candidate window, the present invention merges multiple windows around adjacent SV signals. There are a series of candidate windows on the reference genome, which are composed of SV signals and the regions around them. Each candidate window in the candidate window of the SV signal is represented by Q k =[R s,k ,R e,k ],k=1,2,…,K;

[0129] Among them, R s,k and R e,k They represent the starting position and ending position of the k-th candidate window on the reference genome respectively;

[0130] Each candidate window Q k Corresponding to one or more (greater than or equal to 2) long read sequences, each candidate window Q k The identifiers (each long read sequence in each window has 1 identifier) ​​constitute the set M k ;

[0131] In order to optimize the reliability of the window, the present invention adopts a greedy strategy to merge adjacent windows. The merging conditions include:

[0132] Position proximity: The interval between adjacent windows satisfies Δ k =R s,k+1 -R e,k ≤δ, the adjacent window Q k With Q k+1 Merge into a new window Q′ k =[R s,k ,R e,k+1 ];

[0133] Among them, δ is the distance threshold, indicating that the two windows are close enough; Δ k represents the gap between the starting position of the k+1th window and the ending position of the kth window;

[0134] R s,k+1 represents the starting coordinate of the k+1th window on the reference genome, R e,k Indicates the end coordinate 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. The adjacent window 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 real-time integration and dynamic merging of adjacent windows supported by the same long-read sequence, a series of high-confidence local windows are obtained to integrate scattered SV signals, and 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 the first to sixth embodiments.

[0138] Specific embodiment eight: This embodiment differs from specific embodiments one to seven in that a long read sequence overlap graph is constructed in step three-one; the specific process is as follows:

[0139] Based on the long read sequence set in the remaining window of step 2, the string graph method is used to construct an overlap graph;

[0140] In the initial string diagram (such as Figure 4 ), each non-inclusive sequencing sequence will be split into two vertices - fB and fE, (f refers to the sequencing sequence, fB refers to the starting position of the sequencing sequence, and fE refers to the ending position of the sequencing sequence), corresponding to the positions of the read segments at the 5' end and 3' end respectively. There is no need to explicitly store the overall direction or complementary sequence of the read segment, but only to decide whether to take Watson–Crick complementarity for the read segment according to the vertex type (B or E) during subsequent splicing. 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 kept (see Figure 4 The dual-vertex model can accurately represent the positioning of both ends of the read segment and flexibly handle the splicing relationship of the forward or reverse chain. Figure 4The genomic sequence shown contains two key regions: the yellow-marked regions represent unique sequences in the genome, each of which appears only once in the entire genome, and the red and green regions represent repetitive sequences, where sequences of the same color appear twice and belong to exactly the same repetitive units.

[0141] For the overlap o between any two reads f and g, if the corresponding overlapping 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 reflect the overhanging sequences outside the overlap:

[0142] 1. Left-end extension edge: from vertex f.B to g.B, with the edge label being f[0, o.f.beg], that is, the remaining base sequence before the overlapping region of f;

[0143] 2. Right-end extension edge: from vertex g.E to f.E, with the edge label being comp(g[o.g.end, g.len]), where comp(...) represents Watson–Crick complementarity, that is, taking the complement of the remaining bases after the overlapping region of g.

[0144] When the overlap is in the opposite direction on g (o.g.beg > o.g.end), the same rules apply, except that when selecting the non-overlapping interval, the corresponding section needs to be located first according to the overlap direction. In this way, whether it is a forward or reverse overlap, each edge accurately labels how the two reads are spliced outside the overlap, ensuring that when walking along any path, both the overlapping regions can be automatically aligned and the remaining overhanging sequences and the reads themselves can be pieced together in the correct direction and in accordance with the strand information order;

[0145] The present invention adopts an improved strategy to construct a long-read sequence overlap graph. This method is based on the minimizer technique and can identify the overlapping regions between long-read sequences. In order to avoid missing short SV signals, strict constraints are imposed when constructing the alignment skeleton;

[0146] The overlap graph satisfies the following constraints:

[0147] Sequence overlap length constraint: It is required that the overlapping segments of two long-read sequences within one window are at positions q k+1,start -q k,end ≤Δ max ;

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

[0149] This restriction ensures that adjacent matching blocks are close enough in length sequence, so that continuous matching information can be merged to obtain more accurate information.

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

[0151] Among them, D max It is the preset maximum distance to avoid the introduction of incorrect connections due to the long distance between matching blocks, and to ensure the consistency of the comparison skeleton with the real reference structure; r k+1,start represents the starting coordinates of the k+1th long read sequence overlap fragment on the reference genome; r k,end Indicates the end coordinates of the k-th long read sequence overlap fragment on the reference genome;

[0152] In real-time scenarios, the 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 the first to seventh embodiments.

[0154] Specific embodiment 9: This embodiment differs from any one of specific embodiments 1 to 8 in that in step 32, the overlapping graph is converted into a directed string graph and simplified;

[0155] The specific process is:

[0156] In the second step, the present invention converts the long-read sequence overlap graph into a directed string graph and, through a series of graph simplification operations, generates contigs derived from linear single paths in the graph. This real-time graph construction allows the system to gradually output parts of the graph structure during the assembly process. This stage of operation mainly includes several key steps: first, by removing transitive edges, redundant information in the graph is eliminated, and the structural complexity of the graph is simplified; second, the graph structure is simplified by pruning short branches caused by sequencing errors.

[0157] First, transitive edges are removed to eliminate redundancy and simplify the complexity of the graph. This paper adopts a transitive reduction algorithm to simplify the directed string graph, aiming to simplify the graph representation of long read sequences by accurately pruning redundant connections.

[0158] 1) Delete edges from the overlapping graph obtained in step 3. The specific process is as follows:

[0159] Initialize all nodes (the length sequence includes 1 forward chain and 1 reverse chain, corresponding to 2 nodes) to an unvisited state to ensure that each node is only visited once during the processing;

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

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

[0162] Calculate the longest edge length And set a fault tolerance parameter For each edge e q If there is a path P from u to v (u and v are edges e q The two vertices of It is believed that e q It is a redundant transitive edge, and the edge e q Mark and remove, thereby effectively pruning redundant connections and simplifying the complexity of the string. This process can gradually output the simplified graph structure during real-time graph construction;

[0163] Among them, length(P) represents the length of path P;

[0164] 2) Prune short branch tips caused by sequencing errors from the overlapped graph obtained in 1) to improve graph accuracy. This algorithm is used to prune short tips in the genome assembly graph, which are often generated by sequencing errors. The algorithm removes these potentially misleading structures by identifying and removing short paths connected to a single input and output edge. Short branches caused by sequencing errors are problems with the sequencing itself, which may contain errors, resulting in short branches (short paths that start at the endpoint but cannot continue). The specific process is as follows:

[0165] Assume that the short branch tip caused by sequencing error is a line 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 , then delete the tip to reduce error information and optimize the subsequent decomposition process. The real-time output of the atlas data can help monitor the tip pruning effect;

[0167] If the tip length does not satisfy L tip <L min , then the tip is not deleted;

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

[0169] 3) Obtaining a bubble structure based on the overlapping image obtained in 2), and deleting the overlapping image obtained in 2) based on the bubble structure;

[0170] The specific process is:

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

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

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

[0174] Among them, d(u,v) represents the shortest path length from u to v, D max is 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 and further distinguishes two categories based on the internal path of the bubble:

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

[0177] Among them, 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 the minor haplotype;

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

[0180] 4) Obtain a continuous linear contig.

[0181] The specific process is:

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

[0183] The present invention ensures the accuracy of the splicing path and extracts the global optimal splicing path by the following methods:

[0184] 42) Check whether all nodes on the remaining path in the new overlapping graph obtained in 41) are derived from the same continuous segment of the reference genome, and the relative direction (forward and reverse strands) of two adjacent nodes on the long read sequence is consistent with the relative direction (forward and reverse strands) on the reference genome;

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

[0186] If not satisfied, execute 43);

[0187] 43) In the new overlapping graph obtained in 41), identify the source node with zero in-degree and greater than zero out-degree (the source node is connected to other nodes) and the sink node with in-degree (the number of in-degree edges) greater than zero and out-degree (the number of out-edges) zero (other nodes are connected to the sink node), and use the Dijkstra algorithm to calculate the shortest path from the source node to the sink node, that is, obtain a continuous linear contig; achieve global optimal splicing in extreme repetitive and complex branching environments.

[0188] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0189] Specific embodiment ten: This embodiment differs from one of specific embodiments one to nine in that, in step three-three, the long read sequence alignment information replaces the reference sequence based on the directed string graph simplified in step three-two;

[0190] The specific process is:

[0191] According to the starting and ending positions of the contig in the reference genome, the contig is replaced with the reference sequence;

[0192] Determine the termination condition:

[0193] (1) Mutation detection convergence: Repeat steps 2 to 3. 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) Stability of replacement coverage: When the genome sequence coverage of the replaced region reaches 95% and the replacement fluctuation range (the number of replaced reference sequences accounts for the total number of reference sequences) after repeating steps 2 to 3 twice is less than 2%, it indicates that the key region update has stabilized;

[0195] If any of these conditions are met, the real-time sequencing process is automatically terminated, and the final updated sequence and all alternative sequences are stored and output in a new FASTA file. This dual dynamic termination mechanism ensures efficient analysis and reliable results.

[0196] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0197] Example:

[0198] Experimental setup

[0199] This paper used both real and simulated datasets. The real dataset used a dataset from the widely studied GIAB sample HG002 (NA24385): the PacBio HiFi dataset. The data was then downsampled (5x, 10x, 20x, and 30x) using seqtk, for a total of five real datasets. The simulated datasets used SURVIVOR to randomly simulate 675 and 732 structural variants, respectively, for Arabidopsis thaliana and Drosophila melanogaster. The variant length range was set to 50–10,000 bp, the homozygous proportion of structural variants was set to 0.6, and the coverage was 5x, 10x, 20x, 30x, and 50x, respectively.

[0200] Table 1 Genome information of species and generated dataset information

[0201]

[0202] The method outputs individual genome sequences (FASTA format) and then aligns them to the reference genome using minimap2 (parameters "-ax asm10–cs -r2k"). Sorted and indexed BAM files are generated using samtools. SVs are then called in haploid mode using SVIM-asm, outputting VCF files. Finally, bcftools is used for sorting, bgzip compression, and tabix indexing. The SV detection performance of all five tools was evaluated using Truvari (parameters "-r 1000 -p 0 --dup-to-ins --pickmulti --passonly").

[0203] When evaluating the model's predictive performance, we use metrics including precision, recall, and F1 score. These metrics comprehensively assess 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 true variants among the variant sites predicted by the model. This metric directly reflects the credibility of the variant detection results.

[0206]

[0207] TP indicates that the model correctly predicts a true variant site, while FP indicates that the model incorrectly predicts a normal site. A high accuracy rate means that the variant site reported by the model has a high degree of confidence.

[0208] 2. Recall

[0209] Recall is the proportion of samples that are correctly predicted as positive by the model among all true positive examples. Recall assesses the model's ability to detect true variations.

[0210]

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

[0212] 3. F1 value

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

[0214]

[0215] When both precision and recall are high, the F1 value will also be high; if one of the two is low, the F1 value will be greatly affected, reflecting the imbalance of model performance.

[0216] 4. Experimental results

[0217] In this paper, four most advanced structural variation (SV) detection tools (cuteSV2 v2.1.1, sniffles2v2.2, SVDSS v2.0.0, and SVIM v2.0.0) were selected for comparison. The experimental results are shown in the following table.

[0218] Table 2 Results Precision, Recall and F1 Value

[0219]

[0220]

[0221] The present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A reference genome-based assembly method for single-molecule real-time sequencing, characterized by: The specific process of the method is: Step 1: Build a hash index based on the reference genome; use a phased alignment strategy to build a long-read sequence alignment skeleton; record the long-read sequence alignment skeleton breakpoints as SV signals; the specific process is as follows: Step 1: Build a hash index based on the reference genome; Step 1 and 2: Based on the hash index, a phased alignment strategy is used to build the alignment skeleton of the long read sequence: Phase 1: Build the initial candidate alignment skeleton and output preliminary alignment results in real time; Phase 2: Optimize the alignment backbone of long-read sequences by re-calling duplicate minimizers; Step 13: Record the alignment backbone breakpoints of the long-read sequence as SV signals; Step 2: Cluster the SV signals; the specific process is as follows: Step 21: Construct candidate windows for SV signals; Step 22: Based on the candidate windows 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 31: 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 in step 32, replace the reference sequence with the long read sequence alignment information.

2. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 1, characterized in that: In the step 1, a hash index is constructed based on the reference genome; The specific process is: 1) Given a reference genome, a window size of w and a step size of 1 base; Divide the reference genome into windows starting from the first base of the reference genome to obtain all windows; 2) Use the hash function to calculate the w k-mers in each window. The k-mer with the smallest hash value is defined as the minimizer of the corresponding window. Collect all minimizers of the reference genome; 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 build a hash table H; The key of hash table H is l-mer, and the value of hash table H is the position list of minimizer in linear list L.

3. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 2, characterized in that: In steps 1 and 2, a phased alignment strategy is used to construct an alignment skeleton for long read sequences based on hash indexes: Phase 1: Build the initial candidate alignment skeleton and output preliminary alignment results in real time; Phase 2: Optimize the alignment backbone of long-read sequences by re-calling duplicate minimizers; The specific process is: 1) Construct the initial candidate comparison skeleton and output the preliminary comparison results in real time; the specific process is: 11) For a given long read sequence Q, a window size of w is given and a step size of 1 base is given; Divide the long read sequence Q into windows to obtain all windows; 12) Use the hash function to calculate the w k-mers in each window. The k-mer with the smallest hash value is defined as the minimizer of the corresponding window. Collect the minimizers of all windows of the long read sequence Q. 13) Match the minimizer of each window of the long read sequence Q with the hash index of the reference genome, and locate the candidate alignment position of the minimizer of each window in the long read sequence Q in the reference genome; 14) Let the i-th matching pair after 13) be (q i ,r i ); Among them, q i represents the position of the long read sequence Q corresponding to the i-th matching pair; r i Indicates the position of the reference genome corresponding to the i-th matching pair; 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 )|≤ε′, then the matching pair (q i ,r i ) and (q i+n ,r i+n ) has collinearity, matching pairs (q i ,r i ) and (q i+n ,r i+n ) are merged to form a matching block, matching block MB={(q,r)|q∈[q start ,q end ],r∈[r start ,r end ]}; Among them, ε′ represents the relative position deviation threshold allowed; q i+n represents the position of the long read sequence Q corresponding to the i+1th matching pair, r i+n Indicates the position of the reference genome corresponding to the i+1th matching pair; i represents the sequence number of the current matching pair among all the extracted minimizer matching pairs in 3); n represents the matched pair interval used for collinearity test; 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 end position of q, r start Indicates the starting position of r, r end Indicates the end position of r; 15) Construct a directed acyclic graph and regard each matching block MB as a node in the directed acyclic graph. If the positions of two matching blocks MB on the reference genome are continuous and the order in which the two matching blocks MB appear on the long read sequence Q is consistent with the order in which they appear on the reference genome, then in the i-th matching block MB i To the jth matching block MB j Create a directed edge between them and define the edge weight w ij and penalty weight p ij ; Calculate the score for each MB: in, S(MB j ) represents the jth matching block MB j score; S(MB i ) represents the i-th matching block MB i score; Predecessors represent all the nodes that can be directly connected to the MB i The set of predecessor matching blocks; 16) Select the matching block with the highest score as the last point of the path, and deduce the above formula forward until the first matching block is found to obtain the path, and use the path as the candidate alignment skeleton; Converting the candidate alignment skeleton into an alignment between the candidate sequence and the reference sequence, wherein the alignment information includes the candidate sequence identifier, the position of the candidate sequence in the reference sequence, and alignment skeleton score information; The candidate sequence is the long read sequence currently aligned; 2) Optimize the candidate alignment skeleton; the specific process is: 21) If a minimizer appears more than R times in the reference genome, the minimizer that appears more than R times will not be considered; 22) retaining matching blocks whose alignment skeleton scores exceed a preset threshold value from the candidate alignment skeletons output in real time in 16); 23) For the gap interval between any two adjacent matching blocks, recall all minimizers in the gap interval; 24) Merge the corrected candidate alignment skeletons.

4. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 3, characterized in that: The candidate alignment skeletons after merging the corrections in 24); The specific process is: 241), determine whether the recalled minimizer and the matching block in front of the gap interval maintain collinearity on the reference genome and length sequence; 242), determine whether the recalled minimizer and the matching block after the gap interval maintain collinearity on the reference genome and length sequence; 243), when 2231) and 2232) maintain collinearity at the same time, the recalled minimizer is used as a new matching block; execute 2234); Otherwise, the recalled minimizer is not used as a new matching block; 244), merging the matching blocks retained in the candidate alignment skeleton in 22) with the new matching blocks to generate an ordered matching block list; 245) Determine whether the ordered matching block list generated in 244) needs to be merged to obtain the final matching block list as the merged and corrected comparison skeleton; the specific process is: 2451), sequentially reading adjacent matching block pairs in the ordered matching block list generated in 244); 2452), when the positions of adjacent matching blocks in the reference genome and the positions in the long read sequence overlap, or the length of the gap between adjacent matching blocks is less than the preset merging threshold Δ, the adjacent matching blocks are merged into one matching block; The starting point of the merged matching block in the reference genome takes the minimum starting coordinate of the adjacent matching block in the reference genome, the starting point of the merged long read sequence takes the minimum starting coordinate of the adjacent matching block in the long read sequence, the end point of the merged matching block in the reference genome takes the maximum ending coordinate of the adjacent matching block in the reference genome, and the end point of the merged long read sequence takes the maximum ending coordinate of the adjacent matching block in the long read sequence; 2453) until all adjacent matching blocks are traversed and the merged and corrected comparison skeleton is output.

5. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 4, characterized in that: In the step 1-3, the skeleton breakpoints are recorded as SV signals; The specific process is: Determine whether there are breakpoints in the alignment skeleton after optimization; If there is a breakpoint, the breakpoint is used as a candidate SV signal and the candidate SV signal position information (R s ,R e ,G s ,G e ); while recording the split alignment information; Among them, R s and R e Represent the starting and ending positions of the long read sequence, G s and G e Represent the starting and ending positions on the reference genome, respectively; If there is no breakpoint, record the full-length alignment information.

6. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 5, characterized in that: Constructing a candidate window of the SV signal in step 21; The specific process is: 1) Divide the reference genome into overlapping windows of length Lbp, where the overlapping length of each window is L / 2; As long-read sequencing data continues to be input, the detected SV signals are dynamically assigned to various windows; 2) Select candidate windows containing SV signals through condition 1 or condition 2; Condition 1: For a window Θ in 1), if Then window Θ is regarded as a candidate window; Among them, C sv Indicates the number of SV signals within the window Θ, RD Θ represents the long read sequence coverage depth of window Θ, P Θ It is a custom parameter; Condition 2: For a certain window Θ in 1), as long as the window Θ contains an SV signal, the window Θ is considered as a candidate window.

7. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 6, characterized in that: In the step 22, the candidate windows based on the SV signal are merged with adjacent SV windows; The specific process is: Each candidate window in the candidate window of the SV signal is represented by Q k =[R s,k ,R e,k ],k=1,2,…,K; Among them, R s,k and R e,k They represent the starting position and ending position of the k-th candidate window on the reference genome respectively; Each candidate window Q k Corresponding to one or more long read sequences, each candidate window Q k The identifiers constitute the set M k ; Adjacent windows are merged using a greedy strategy. The merging conditions include: Position proximity: The interval between adjacent windows satisfies Δ k =R s,k+1 -R e,k ≤δ, the adjacent window Q k With Q k+1 Merge into a new window Q′ k =[R s,k ,R e,k+1 ]; Among them, δ is the distance threshold; Δ k represents the gap between the starting position of the k+1th window and the ending position of the kth window; R s,k+1 represents the starting coordinate of the k+1th window on the reference genome, R e,k Indicates the end coordinate of the k-th window on the reference genome; Long read sequence connectivity: adjacent windows are connected by the same long read sequence, i.e. The adjacent window 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 .

8. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 7, characterized in that: In the step 31, a long read sequence overlap graph is constructed; The specific process is: Based on the long read sequence set in the remaining window of step 2, the string graph method is used to construct an overlap graph; The overlapping graph satisfies the following constraints: Sequence overlap length constraint: requires that the overlapping fragments of two long read sequences within a window be q on the long read sequence coordinates. k+1,start -q k,end ≤Δ max ; Among them, Δ m ax is the preset maximum allowed interval; q k+1,start Indicates the starting position of the overlapping fragment of the k+1th long read sequence; q k,end Indicates the end position of the overlapping fragment of the k-th long read sequence; Reference genome position overlap constraint: It requires that the coordinate intervals of the overlapping fragments of two long read sequences within a window on the reference genome must overlap, and the distance between the coordinates of the overlapping fragments of the two long read sequences on the reference genome must satisfy |r k+1,start -r k,end |≤D max ; Among them, D max is the preset maximum distance; r k+1,start represents the starting coordinates of the k+1th long read sequence overlap fragment on the reference genome; r k,end Indicates the end coordinates of the k-th long read sequence overlap fragment on the reference genome.

9. The reference genome-based assembly method for single-molecule real-time sequencing according to claim 8, characterized in that: In the step 32, the overlapping graph is converted into a directed string graph and simplified; The specific process is: 1) Delete edges from the overlapping graph obtained in step 3. The specific process is: Initialize all nodes to an unvisited state; Let the outgoing edge set of a node u in the overlapping graph be E out (u)={e1,e2,…,e q ,…,e n }; Among them, each edge e q The length of l(e q ); Calculate the longest edge length And set a fault tolerance parameter For each edge e q If there exists a path P from u to v that satisfies It is believed that e q It is a redundant transitive edge, and the edge e q Mark and remove; Among them, length(P) represents the length of path P; 2) Prune the short branch tips caused by sequencing errors from the overlap graph obtained in 1); The specific process is: Assume that the short branch tip caused by sequencing error is a line from endpoint v1→v2→…v q →…→v n The edge sequence {e1,e2,…,e q ,…,e n }, the total length of the tip is If the tip length satisfies L tip <L min , then delete the tip; If the tip length does not satisfy L tip <L min , then the tip is not deleted; Among them, L min Indicates the preset threshold; v1 indicates the starting node of the tip, v q Represents the first intermediate node in the tip path, v n Indicates the terminal node of the tip; 3) Obtaining a bubble structure based on the overlapping image obtained in 2), and deleting the overlapping image obtained in 2) based on the bubble structure; The specific process is: 31) The assembly graph is a directed graph G = (V, E), where V is the node set and E is the edge set; 32) Use depth-first search to traverse the graph G. For each node v∈V, if the out-degree of node v is greater than 1, then v is considered to be a potential bubble entrance; otherwise, v is not considered to be a potential bubble entrance. 33) For two candidate nodes u, v∈V, if there are at least two mutually disjoint paths P1 and P2 from u to v, and satisfy d(u,v)=min{length(P1),length(P2)} <D max , then the paths P1 and P2 structures are defined as bubble structures B(u,v); otherwise, the paths P1 and P2 structures are not bubble structures B(u,v); Among them, d(u,v) represents the shortest path length from u to v; D max is the preset maximum detection distance; length(P1) represents the length of path P1; length(P2) represents the length of path P2; 34) After the bubble structure B(u,v) is detected, the internal path of the bubble is further distinguished into two categories: The main bubble path is denoted as P main =argmaxS(P); Where S(P) represents the number of SV signals contained in all paths P in the bubble structure B(u,v); The alternative bubble path is denoted as P alter ,S(P alter ) <S(P main ); Among them, S(P main ) represents the main bubble path P main The number of SV signals contained in ; 4) Obtain a continuous linear contig; The specific process is: 41) Generate an alternative sequence based on the alternative bubble path, delete the alternative sequence sequence in the overlapping graph obtained in 3) to obtain a new overlapping graph; 42) Check whether all nodes on the remaining paths in the new overlapping graph obtained in 41) are derived from the same continuous segment of the reference genome, and whether the relative directions of two adjacent nodes on the long read sequence are consistent with the relative directions on the reference genome; If satisfied, a continuous linear contig is obtained; If not satisfied, 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, and use the Dijkstra algorithm to calculate the shortest path from the source node to the sink node, that is, obtain a continuous linear contig.

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

Images