Genome stacks, genome architecture, genome sequence assembly methods and systems

Abstract

This invention provides a method and system for assembling a genome stack, a method and system for assembling a genome architecture, and a method and system for assembling a genome sequence, belonging to the field of bioinformatics. This invention optimizes a global loss function using robust regression techniques to determine the stacking relationships between sequencing sequences, eliminating interference from false positive alignments, thereby obtaining a more accurate assembled stack set and avoiding copy number loss or misassembly in the assembled genome. The estimation of stack arrangement based on regression matrix calculations has a higher tolerance for outliers; unless the number of ambiguous sample information exceeds the number of true samples, a "collapse" will not occur. By using resampling techniques to independently sample, assemble, improve, and evaluate sequencing data multiple times, and then integrating them into the final assembled genome, the uncertainty of the assembly results caused by noise in the sequencing data can be reduced, and the sensitivity of the assembly results to the selection of assembly parameters can be decreased.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to a genome array assembly method, a genome architecture assembly method, a genome sequence assembly method, a genome array assembly system, a genome architecture assembly system, and a genome sequence assembly system. Background Technology

[0002] Genome sequencing is a crucial technique for conducting molecular biology research. Researchers can obtain the genomic base sequence of an organism through genome sequencing. This sequence serves as a genetic template, providing an important reference for studying gene transcription, regulation, and modification, thereby helping to explain the molecular mechanisms behind life phenomena. By comparing the genomes of different species or individuals, researchers can discover differences at the genomic level, providing guidance for revealing evolutionary history, improving crop breeding, diagnosing genetic diseases, and optimizing drug treatments.

[0003] Currently, the most widely used sequencing technologies are second-generation (NG) and third-generation (NG) sequencing technologies. NNG sequencing offers high base identification quality and is inexpensive, typically producing sequences (i.e., the number of base pairs) between 100-300 bp. A key feature is its ability to sequence very long fragments from both ends, obtaining the base sequences at both ends of the fragment—a paired-end sequence. Third-generation (NG) sequencing technologies can measure significantly longer sequences (typically between 10-100 kbp), but this is accompanied by increased sequencing costs and a higher base identification error rate.

[0004] Assembling sequenced sequences into a genome is a fundamental problem in computational biology. Because the sequence length that sequencers can measure is much smaller than the length of the genome, it is necessary to assemble the sequences after sequencing to infer their relative positions and thus reconstruct the sequenced genome. A key challenge in genome assembly is that the genome contains many highly similar or repetitive regions, the presence of which greatly increases the uncertainty in inferring the relative positions of sequenced sequences.

[0005] Existing genome assembly methods are mainly divided into two categories in principle. One category is based on the De Bruijn diagram. The main operation of this method is as follows: for each sequencing sequence, a subsequence of a specific length (usually called a k-mer, where k represents the length of the subsequence) is cut off every other base; a De Bruijn diagram is constructed using all the cut subsequences; after certain error correction operations, paths are found on the diagram, and each path is inferred to be a fragment on the genome being analyzed. The other category is based on the stacking relationship of sequencing sequences. This method compares every two sequences, and then infers the stacking relationship of the sequences based on the alignment results, thereby obtaining an assembled stack set, where each stack set corresponds to a fragment on the genome being analyzed.

[0006] After obtaining several stacks (or fragments), they need to be oriented and arranged to assemble into a structure. The general method is to compare the sequencing sequence with all the congruent sequences of the stacks, and use the mapping and library length information of the paired-end sequencing sequence, or the coordinate information of the single-molecule long-read sequencing sequence mapped to multiple congruent sequences of the stacks, to determine the order, orientation relationship and distance range between different stacks, thereby obtaining the assembled genome structure.

[0007] The De Bruijn diagram-based method works well for assembling genomes with low repetition, but is not ideal for genomes with high repetition because the length of the cut subsequences is significantly shorter than the sequencing sequence, thus reducing specificity and causing errors when finding paths on the diagram.

[0008] Methods based on sequence alignment use sequencing sequences as units, rather than k-mers, thus preserving more positional information across the genome and facilitating the reconstruction of repetitive sequences. However, current methods primarily rely on a greedy strategy, sequentially selecting locally optimal alignments to determine sequence alignment. These methods are susceptible to false positives, leading to incorrect alignments and resulting in copy number loss or misassembly in the assembled genome. While using stricter sequence alignment criteria to reduce false positives can increase false negatives, it can also decrease the continuity of the assembled genome.

[0009] Because the length of the library in paired-end sequencing can sometimes deviate significantly, and because repetitive sequences and heterozygous regions on the genome can cause ambiguity in the arrangement of stacks, existing genome architecture methods struggle to achieve a genome architecture that balances length and accuracy when determining the orientation and relative position between stacks.

[0010] Due to deviations in sample preparation and base identification, the sequencing sequences may contain errors inconsistent with the genome being analyzed, leading to inaccurate assembly results. For the same genome assembly method and the same genome being analyzed, different input sequencing datasets often produce different assembly results. Existing genome assembly systems typically use all sequencing sequences as input, failing to eliminate the uncertainty introduced by sequencing data into the assembly results. Furthermore, the assembly results often depend on the selection of assembly parameters, and parameter tuning is usually based on experience. Summary of the Invention

[0011] The purpose of this invention is to provide a method and system for assembling genome stacks to eliminate the misleading effect of false positive alignments in determining the stacking relationship of sequencing sequences, thereby increasing the accuracy and continuity of genome assembly.

[0012] To achieve the above objectives, this invention provides a genome stacking array assembly method. The method first randomly extracts a first sequencing sequence set at a certain depth from the sample sequencing sequences; then, it performs pairwise alignments of the sequencing sequences in the first sequencing sequence set using a predetermined sequence alignment algorithm to obtain pairwise alignment information; based on the pairwise alignment information, all alignments satisfying the stacking condition are collectively represented as a first linear regression model; the first linear regression model is solved using an iterative reweighted least squares algorithm to obtain the solution with the minimum global loss function and a robust estimate of the coordinates of the sequencing sequences in the first sequencing sequence set on the genome being tested; the sequencing sequences are divided into multiple connected components using an undirected graph constructed based on the solution with the minimum global loss function, and the sequencing sequences are arranged according to the robust estimates of the coordinates of the sequencing sequences within each connected component to generate a stacking array; the stacking arrays are merged based on the pairwise alignment information of the sequencing sequences to obtain a preliminary assembled stacking array set; the consensus sequence of each stacking array in the preliminary assembled stacking array set is obtained, and the sequencing sequence alignments are mapped to the consensus sequences to obtain the mapping result from the sequencing sequences to the consensus sequences, which is the new assembled stacking array set.

[0013] The present invention also provides a genome array assembly system, which uses the aforementioned genome array assembly method to assemble genome arrays.

[0014] Another objective of this invention is to provide a method and system for assembling genome architectures to eliminate the interference of ambiguous information in the stacked array arrangement on the architecture assembly, thereby balancing the length and accuracy of the assembled architecture.

[0015] To achieve the above objectives, the present invention provides a genome architecture assembly method, which first randomly extracts a second sequencing sequence set at a certain depth from the sample sequencing sequence; then assembles a genome stack based on the second sequencing sequence set to obtain the assembled stack set. The second sequencing sequence set and the assembled stack set The sequence set is aligned with the target genome, and single-mapped sequencing sequence stacks are selected based on the alignment results to obtain uniquely mapped sequencing sequence information and the selected sequence set. An undirected weighted stacking matrix is ​​constructed based on the uniquely mapped sequencing sequence information and the selected sequence set, and the undirected weighted stacking matrix is ​​divided into multiple connected components. A combinatorial optimization model is used to determine the orientation of the stacking matrix corresponding to each vertex in each connected component. A regression model is established to estimate the arrangement coordinates of the stacking matrix on the genome being tested. All stacking matrices in the same connected component are sorted according to the orientation of the stacking matrix corresponding to each vertex and the arrangement coordinates to obtain the quasi-architecture of that connected component. An architecture set is obtained based on the quasi-architectures of multiple connected components. The arrangement coordinates of the stacking matrix in each architecture in the architecture set are sorted to eliminate overlapping or redundant stacking matrices, resulting in the assembled genome architecture.

[0016] The present invention also provides a genome architecture assembly system, which uses the aforementioned genome architecture assembly method to assemble genome architectures.

[0017] Another objective of this invention is to provide a genome sequence assembly method and system to reduce the uncertainty of assembly results caused by noise in sequencing data and to reduce the sensitivity of assembly results to the selection of assembly parameters.

[0018] To achieve the above objectives, the present invention provides a genome sequence assembly method. This method first randomly extracts a third sequencing sequence set at a certain depth and a fourth sequencing sequence set at another depth from a sample sequencing sequence. The third sequencing sequence set is then processed using the genome stack assembly method to obtain an assembled stack set. The assembled stack set is then processed using the genome architecture assembly method to obtain an assembled genome architecture. The assembled genome architecture is updated and improved using the fourth sequencing sequence set, and the improved assembled genome is quality-assessed. If the quality meets the standard, the improved assembled genome is stored as a candidate assembled genome. The aforementioned steps are repeated until a predetermined number of candidate assembled genomes are obtained. All candidate assembled genomes are then subjected to multiple sequence alignment, and the most frequent bases in each aligned column are linked together to obtain the final assembled genome.

[0019] The present invention also provides a genome sequence assembly system, which uses the aforementioned genome sequence assembly method to assemble genome sequences.

[0020] The above technical solutions can eliminate the misleading effect of false positive alignments in determining the stacking relationship of sequencing sequences, thereby increasing the accuracy and continuity of genome assembly. At the same time, they can eliminate the interference of ambiguous information in the stacking arrangement on the architecture assembly, thus balancing the length and accuracy of the assembled architecture. They can also reduce the uncertainty of the assembly results caused by noise in the sequencing data and reduce the sensitivity of the assembly results to the selection of assembly parameters.

[0021] Other features and advantages of the embodiments of the present invention will be described in detail in the following detailed description section. Attached Figure Description

[0022] The accompanying drawings are provided to further illustrate embodiments of the present invention and form part of the specification. They are used together with the following detailed description to explain the embodiments of the present invention, but do not constitute a limitation thereof. In the drawings:

[0023] Figure 1 This is a schematic diagram illustrating the establishment of a linear regression model using pairwise alignment information of sequencing sequences according to an embodiment of the present invention;

[0024] Figure 2 This is a schematic diagram illustrating the necessary merging of connected components derived from regression results based on pairwise alignment information of sequencing sequences according to an embodiment of the present invention;

[0025] Figure 3 This is a schematic diagram illustrating the remapping of sequencing sequences to congruent sequences according to an embodiment of the present invention;

[0026] Figure 4 This is a general flowchart of a method for assembling a genome architecture according to an embodiment of the present invention;

[0027] Figure 5 This is a schematic diagram illustrating the extraction of pairwise directional relationships and start-point distance information between identical sequences from paired-end sequencing sequence mapping data according to an embodiment of the present invention.

[0028] Figure 6 This is a schematic diagram of extracting the directional relationship and start-point distance information between pairs of congruent sequences from single-molecule long-read sequencing sequence mapping data according to an embodiment of the present invention.

[0029] Figure 7 This is an algorithm flowchart for optimizing the direction assignment of each vertex on a connected component of a stacked matrix using a combinatorial optimization model, according to an embodiment of the present invention.

[0030] Figures 8a-8b This is a schematic diagram illustrating eight scenarios for calculating the directional relationships and distance information between pairs of quasi-archites according to an embodiment of the present invention;

[0031] Figure 9This is a schematic diagram illustrating the final architectural arrangement based on regression-estimated coordinates, connectivity information, and local comparison of overlapping regions according to an embodiment of the present invention.

[0032] Figure 10 This is a block diagram of a resampling-based genome sequence assembly system according to an embodiment of the present invention. Detailed Implementation

[0033] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0034] To facilitate understanding of the genome sequence assembly method in the embodiments of the present invention, the relevant technical terms in the embodiments of the present invention are defined and explained as follows:

[0035] Genome: A genetic sequence containing all the genetic information of an organism. It is composed of four bases and is usually represented as a sequence containing four characters: A, C, G, and T, each of which represents a base.

[0036] Genome to be tested: The genome whose base sequence needs to be determined.

[0037] Sequencing sequence: The measurement result of a subsequence of an organism's genome output by a sequencing platform. It consists of several characters representing the four bases (A, C, G, T), representing a subsequence of the genome. One end of the sequencing sequence is the 5' end, and the other end is the 3' end; usually, in the file provided by the sequencing platform, the left end of the sequencing sequence is the 5' end, and the right end is the 3' end.

[0038] bp: Abbreviation for basepair. 100bp means a sequence of 100 bases. 1kbp means a sequence of 1000 bases.

[0039] Depth: The ratio of the total number of bases in the sequenced sequence to the size of the genome being sequenced, which is the average number of times a single base on the genome is read. For example, if the sequencing depth of a sample is 30X, it means that each single base on the genome of that sample was read an average of 30 times.

[0040] Second-generation sequencing technology, also known as next-generation sequencing technology, is characterized by its ability to generate a large number of paired-end sequencing data in parallel, with each sequencing sequence typically ranging from 100 to 300 bp in length.

[0041] Paired-end sequencing sequence: In the second-generation sequencing technology, the sequence to be tested is cut into several fragments, the length of which can be small or large. Then, sequencing is performed at both ends of the fragments, resulting in a sequencing sequence at one end and a sequencing sequence at the other end. If the fragment is too long, the middle part is generally not detectable. The two sequencing sequences obtained correspond to the same fragment and are called paired-end sequencing sequences.

[0042] Fragment library and library length: In the same paired-end sequencing sequence file provided by the sequencing platform, all paired-end sequencing sequences belong to fragments with similar lengths, and the sequencing platform will estimate an average value; in the following text, each paired-end sequencing sequence file provided by the sequencing platform is called a fragment library, and its estimated average fragment length is called the library length.

[0043] Short library, long library: In paired-end sequencing, if the sequence fragment to be sequenced is first circularized and then broken up for sequencing, the resulting fragment library is called a long library (English name: mate-paired library), which is usually 2-10kbp in length, and the sequencing direction is from the inside of the fragment to both ends; if the sequence fragment to be sequenced is not circularized, the resulting fragment library is called a short library (English name: paired-end library), which is usually 200-500bp in length, and the sequencing direction is from both ends of the fragment to the inside.

[0044] Third-generation sequencing technology refers to single-molecule sequencing technology, also known as Third Generation Sequencing, Single Molecule Sequencing, or Long-read Sequencing. Currently, it mainly includes two technical directions: single-molecule fluorescence sequencing and nanopore sequencing. Its characteristic is that it can generate long-read sequencing data, with each sequencing sequence being approximately 10-100 kbp in length, known as single-molecule long-read sequencing sequences. However, the sequencing error rate is relatively high, generally around 10%.

[0045] Quality score: During sequencing, each base in the sequence has the potential to be misread. The quality score is a measure of this potential. Each base in the sequence corresponds to a quality score; the higher the quality score, the lower the potential for misreading.

[0046] Insertion: refers to the insertion of an extra base sequence between two adjacent bases in the sequencing sequence relative to the genome being tested.

[0047] Cleavage: refers to the loss of one or more base sequences relative to the genome being sequenced.

[0048] INDEL: An insertion or deletion operation.

[0049] Base substitution: After the sequencing sequence is mapped onto the reference genome, some bases are different from their corresponding bases on the reference genome.

[0050] Sequencing sequence mapping: For a given sequencing sequence, a subsequence of roughly the same length is found on a reference genome. This subsequence is either identical to the sequencing sequence or differs from it, but within predefined criteria. The difference is typically measured by the number of base substitutions or the length of the INDEL.

[0051] Mapping Results: For a sequencing sequence, if a subsequence exists on the reference genome and the difference between the subsequence and the sequencing sequence is within the set standard, then the sequencing sequence is successfully mapped. The mapping results include: (1) the coordinates of the rightmost base of the subsequence on the reference genome, which serves as the mapping coordinates of the sequencing sequence on the reference genome; (2) the mapping direction of the sequencing sequence, which can be mapped in the direction of 5' to 3' or in the direction of 3' to 5'; (3) the alignment information between the sequencing sequence and the subsequence on the reference genome, including base substitutions and INDEL conditions between them. When there are multiple subsequences on the reference genome that differ from the sequencing sequence within the set standard, the sequencing sequence has multiple mapping results, each mapping result corresponding to a subsequence on the reference genome.

[0052] Genome assembly: The genome being analyzed is reconstructed using sequencing sequences, and the resulting genome is called an assembled genome.

[0053] Sequencing sequence contigs (or simply contigs): The English term is contig or layout; it represents a collection of sequencing sequences where the stacking relationship or relative distance between any two sequences is fixed. A sequencing sequence contig can be represented by a matrix containing five characters: A, C, G, T, and Ф. An element of the matrix with the digit Ф indicates that there is no base at that position. Each column of the matrix corresponds to a position in the contig and contains at least one non-Ф character. The number of columns in the matrix is ​​the length of the contig. Each row of the matrix represents a sequencing sequence, and the column index of the last non-Ф character in each row is the coordinate of the sequence in the contig within the contig.

[0054] Stacked array set: A set composed of stacked arrays.

[0055] Consensus (or lattice consensus sequence): The English term is consensus; for each column of the lattice, a base is inferred from all the base information contained therein, and the inferred bases are concatenated to obtain the sequence; a feasible method for inferring bases is to take the base with the highest frequency among all the bases contained in each column of the lattice.

[0056] Consistent sequence set (overlapping consistent sequence set): A set consisting of consistent sequences.

[0057] Stacked matrix graph: A graph constructed using stacked matrix consistent sequences as vertices and edges connecting the vertices based on the mapping results from sequencing sequences to stacked matrix consistent sequences.

[0058] Architecture (or assembled genome architecture): The English term is scaffold; it consists of a series of stacked arrays arranged in a certain order, where the orientation and relative positional relationships of all stacked arrays are known, serving as a framework for assembling the genome.

[0059] Architecture set: A collection of architectures.

[0060] Gap: The English term is gap; the unassembled blank area between adjacent stacks in the architecture whose order is determined. In the assembly result, it is generally filled with N, where the number of N represents the estimate of the gap length.

[0061] This invention uses paired-end sequencing sequences or single-molecule long-read sequencing sequences as inputs, wherein the paired-end sequencing sequences and single-molecule long-read sequencing sequences are derived from the genome being tested, and the paired-end sequencing sequences can contain multiple fragment libraries with different library lengths.

[0062] This invention provides a method for assembling a genome stack, the method comprising:

[0063] S01: Extract a first sequencing sequence set at a certain depth from the input sample sequencing sequence, with a depth of d1. The sequences in the first sequencing sequence set are denoted as follows:

[0064] S02: The sequencing sequences in the first sequencing sequence set are compared pairwise using a sequence alignment algorithm to obtain pairwise alignment information. Before alignment, the parameters used for alignment and the criteria for judging whether the alignment score is significant can be set according to the length and quality distribution of the sequencing sequences, as well as the requirements for alignment sensitivity and specificity.

[0065] S03: Based on the pairwise alignment information of the sequencing sequences, all alignments that meet the overlap condition are collectively represented as the first linear regression model. Specifically:

[0066] S031: The overlap condition is set as follows: the alignment score reaches the significance standard, and the length of the unaligned subsequences at both ends of the overlap region is less than a given threshold.

[0067] S032: Establish two coordinate axes along the two single strands of the DNA sample from the 5' end to the 3' end; any sequencing sequence R from the first sequencing sequence set. i Its coordinates on the single strand of DNA are defined as R. iThe coordinate of the rightmost base is denoted as β. i Similarly, R can be defined. i The reverse complementary sequence (denoted as Coordinates on complementary single chains R is the parameter to be estimated; for any two sequencing sequences that satisfy the overlap condition, R i and R j Assuming the sequencing sequence R i Align the rightmost base to the sequencing sequence R j Internally, R j The right overhang contains y i,j One base, thus obtaining an observation sample: y i,j =β j -β i +ε i,j ; where ε i,j β represents the observation error of the observed sample. j It is R j The coordinates on the single strand of DNA it is located at;

[0068] S032: Assuming the second sequencing sequence R j The leftmost base is aligned to the first sequencing sequence R. i Internally, R i The left overhang contains Each base pair yields a dual observation sample: in, This represents the observation error of the dual observation sample. It is R j The coordinates of the inverse complementary sequence on the complementary single chain;

[0069] S033: Based on the pairwise alignment information of the sequencing sequences, all aligned observation samples and their pairs that meet the overlap condition are integrated into a matrix equation: Y = Xβ + ε; where each row represents an observation sample, and the vector... For the parameters to be estimated, X is a large sparse matrix. Each row of X has only two non-zero elements, 1 and -1. X is stored in sparse matrix format.

[0070] It should be noted that, in this embodiment, the observation error ε i,j Specifically, it refers to R j The right overhang portion is calculated by subtracting the number of deleted bases from the number of bases inserted relative to the tested genome; similarly, This represents the observation error of the dual observation sample, specifically R0. i The left overhang represents the number of bases inserted relative to the tested genome minus the number of bases deleted; ε in the matrix equation represents the vector (ε1,...,ε) m )T , where m represents the total number of observation samples and dual observation samples.

[0071] In some other embodiments, the coordinates of the sequencing sequence can be defined as the coordinates of the leftmost base of the sequencing sequence, then the observed sample is represented as: y i,j =β j -β i +ε i,j ; where y i,j R represents i The number of bases in the left overhang portion, ε i,j R represents i The left overhang portion is the number of bases inserted relative to the number of bases deleted from the number of bases inserted into the tested genome.

[0072] Figure 1 This is a schematic diagram illustrating the establishment of a linear regression model based on sequencing sequence alignment information, as shown below. Figure 1 As shown, each pair of sequencing sequences that satisfy the overlap condition provides one observation sample for the model, while the reverse complementary sequence of these two sequences provides another dual observation sample for the model.

[0073] S04: Solve the first linear regression model using an iterative reweighted least squares algorithm to obtain the solution with the minimum global loss function and a robust estimate of the coordinates of the sequencing sequences in the first sequencing sequence set on the genome being tested, including:

[0074] S041: Solving sparse linear equations about β using accelerated algorithms (e.g., LGMRES algorithm) for large sparse linear equation systems: X T X·β=X T Y, the solution obtained is The initial value;

[0075] S042: In each iteration, based on the current... Calculate the residual vector as well as Maximum value of the absolute value of each component

[0076] S043: If Then calculate the weight matrix. Where α is a preset convergence threshold; the higher the base identification error rate of the sequencing sequence, the larger the value of α, and vice versa; the weighting function is...

[0077] S044: Based on the calculated W, the smaller of the weight values ​​of the two samples that are each other is assigned to these two samples. That is, assuming that the weights of the two samples that are each other are a1 and a2, where a1 > a2, then a2 is used as the weight of these two samples.

[0078] S045: Solving sparse linear equations about β using accelerated algorithms (e.g., LGMRES algorithm): X T WX·β=X T WY;

[0079] S046: Update using the obtained solution Delete the rows in matrices X and Y with weights of 0, and return to the next iteration;

[0080] S047: Until Stop the iteration and record the matrix X from the last iteration;

[0081] S048: For the final round Output after rounding each component This is a robust estimate of the coordinates of the sequencing sequences in the first sequencing sequence set on the genome being tested.

[0082] The function ρ can be any function that simultaneously satisfies the following three conditions, and it has achieved good results in practical applications:

[0083] (1) ρ(x)≥0 holds for all χ, and ρ(x) takes the minimum value at χ=0;

[0084] (2) ρ(x) is symmetric about the Y-axis;

[0085] (3) When χ≥0, ρ(x) is monotonically non-decreasing with respect to χ, and ρ(x) will not become too large when χ increases.

[0086] In some embodiments, the function ρ employs a dual-weight loss function ρ B :

[0087] The corresponding weighting function is:

[0088]

[0089] In this invention, parameter c should be about A monotonically increasing function, and

[0090] In some embodiments, the function that yields ideal results after testing is the following step function:

[0091]

[0092] S05: Based on the undirected graph constructed using the solution that minimizes the global loss function, the sequencing sequence is divided into multiple connected components. The sequencing sequences are then arranged according to robust estimates of the coordinates within each connected component to generate a stacked matrix, including:

[0093] S051: Using the sequencing sequences in the first sequencing sequence set and their corresponding reverse complementary sequences as vertices, a total of 2n1 vertices are contained, each representing a sequencing sequence. Two sequencing sequences are determined based on the two non-zero elements in each row of the matrix X in the last iteration. An edge is established between the vertices corresponding to these two sequencing sequences to construct an undirected graph.

[0094] S052: Traverse the undirected graph and divide it into multiple connected components, obtaining 2k connected components, denoted as C1, ..., C2k. k , Where C i and The included sequencing sequences are one-to-one correspondences and are inverse complementary sequences; if a connected component simultaneously includes a certain R... i and Then remove this connected component;

[0095] S053: Arrange the sequencing sequences in each connected component according to a robust estimate of their coordinates on the genome being sequenced, resulting in a stacked matrix: L1, ..., L k , Stacked array L i and These correspond to two subsequences on the DNA double helix that are inversely complementary.

[0096] It should be noted that, generally, sequencing sequences are arranged from left to right in ascending order. However, in some other embodiments, sequencing sequences can be arranged from right to left in descending order. Taking the connected component C1 as an example, the sequencing sequences within it are arranged sequentially in descending order of their coordinate estimates. For any two adjacent sequencing sequences, the number of bases separating their right ends is equal to the difference in their coordinate estimates. The resulting arrangement is the stacked array L1.

[0097] S06: Based on the pairwise alignment information of the sequencing sequences, the stacked arrays are merged to obtain a preliminary assembled stacked array set, including:

[0098] S061: The set of stacked arrays includes 2k vertices, representing stacked arrays L1, ..., L2 respectively. k ,

[0099] S062: For any two superimposed matrices, without loss of generality, Li and L j If L i The 3' end s sequencing sequences and L j If there are more than t alignments satisfying the overlap condition among s sequencing sequences at the 5' end, then an L-band is established. i The corresponding vertex points to L j Construct a directed graph using the directed edges corresponding to the vertices; where the weight of the edge is set to the number of logs that satisfy the stacking condition, and the values ​​of s and t are positively correlated with the sampling depth d1;

[0100] S063: From the directed graph, the stack L containing the most sequencing sequences x Begin by determining the starting point from the stacked matrix L. x The stacked matrix L pointed to by the edge with the largest weight at the starting point y ;

[0101] S064: From L x A sequencing sequence R is selected from the 3' end. i , and from L y A sequencing sequence R is selected from the 5' end. j R i and R j The comparison between them satisfies the overlap condition;

[0102] S065: Assume R i The rightmost base is aligned to R j Internally, R j If the right overhang contains Δ bases, then L y Robust estimates of the coordinates of all sequencing sequences are added Assume L y Any sequencing sequence R in u Its coordinates are estimated by Become

[0103] S067: L x and L y The sequencing sequences in the array are arranged and merged into a stacked array according to robust estimates of their coordinates from smallest to largest, while L... x and Remove the corresponding vertex in the directed graph;

[0104] S068: From L y Repeat steps S063-S067 until no more stacks can be merged. Then, starting from the stack with the most sequencing sequences among the remaining vertices in the directed graph, repeat steps S063-S067 until all vertices in the directed graph have been traversed. The resulting merged stacks are taken as the initial assembled stack set, denoted as D1,...,D m .

[0105] Figure 2 This is a schematic diagram illustrating the merging of stacked arrays based on sequencing sequence alignment information, as shown below. Figure 2 As shown, the three stacked matrices are merged into a single stacked matrix by the connection information provided by the alignment between the 5' and 3' end sequences.

[0106] S07: Obtain the consistent sequences of each array in the preliminary assembled array set, map the sequencing sequences to the consistent sequences, and obtain the mapping result from the sequencing sequences to the consistent sequences. The mapping result is a new assembled array set, including:

[0107] S071: For any one of the stacked arrays D in the preliminary assembled stacked array set i Take the base with the highest frequency in each column to form a stacked array D. i The corresponding consistent sequence is denoted as E. i ;

[0108] S072: For the stacked array D i Any sequencing sequence R in j The leftmost base is at D i The coordinates in the diagram are a1, R j The rightmost base is at D i The coordinates in E are marked as a2. i Extracting the subsequence from the a1-τth base to the a2+τth base and comparing it with R j Perform local sequence alignment;

[0109] S073: The portion of the local sequence that matches the alignment is taken as R. j To E i The mapping result, D i All sequencing sequences to E i The mapping results constitute a new superposition matrix. For the updated assembled array set.

[0110] Figure 3 This is a schematic diagram illustrating the remapping of sequencing sequences to consistent sequences, as shown below. Figure 3 As shown, for any sequencing sequence in the stacked array, its corresponding position in the conformance sequence is found, and the subsequence of the conformance sequence at the position, along with the base sequences of a predetermined length on both sides, is extracted and compared with the sequencing sequence.

[0111] A second aspect of the present invention provides a genome array assembly system, which performs genome array assembly using the aforementioned genome array assembly method.

[0112] The genome stacking assembly method provided in this embodiment uses robust regression technology to optimize a global loss function to determine the stacking relationship between sequencing sequences. This can eliminate the interference of false positive alignments, thereby obtaining a more accurate assembled stacking set and avoiding copy number loss or incorrect splicing on the assembled genome. The optimization algorithm is based on solving a large sparse linear equation system, with sufficient theoretical research to support the economical use of computation time and storage space, and can be accelerated by GPU or distributed computing. The estimation of the stacking arrangement is based on regression matrix calculation rather than graph calculation, which has a higher tolerance for outliers. Unless the number of ambiguous sample information is greater than the number of real samples, there will be no "collapse" situation: that is, the estimated arrangement value differs greatly from the actual arrangement.

[0113] This invention also provides a method for assembling a genome architecture, such as... Figure 4 As shown, the method includes:

[0114] S11: Extract a second sequencing sequence set at a certain depth from the input sample sequencing sequence, with a depth of d1. The sequences in the second sequencing sequence set are denoted as follows:

[0115] S12: Assemble the genome stack based on the second sequencing sequence set to obtain the assembled stack set.

[0116] S13: Transfer the second sequencing sequence set R1, R2, ..., R... n1 With assembled stacked array The sequence set is compared with the target sequence set. Based on the comparison results, single-mapped sequencing sequence arrays are selected. The single-mapped sequencing sequences and their mapping arrays are then filtered according to set conditions to obtain uniquely mapped sequencing sequence information and the selected target sequence set. Specifically,

[0117] If the input sequencing sequence is a paired-end sequencing sequence, then:

[0118] S1301: Select the mapping results from the alignment results, mapping both ends of the paired-end sequences to the same congruent sequence, to estimate the library length of the paired-end sequencing fragment library. The mapping results include the coordinates of the starting point of each end sequence alignment to the congruent sequence, denoted as s1,t1 and s2,t2 respectively. The observed value of the library length is denoted as: in = max(s1,t1,s2,t2) - min(s1,t1,s2,t2). For each fragment library, a predetermined number of observations are taken, and the median M of the observations is taken as the observed library length, and its standard deviation σ is calculated. In one embodiment of the present invention, 30,000 observations are taken for each fragment library, and the median of these 30,000 observations is taken as the observed library length.

[0119] S1302: The mapping coordinates of the paired-end sequencing sequences mapped to two different stacked consistent sequences are still denoted as s1,t1 and s2,t2, respectively.

[0120] S1303: Determine whether the direction of the paired-end sequencing sequence and its mapped congruent sequence is consistent based on the mapped coordinates: Due to the principles of next-generation sequencing, the sequencing direction of a short library is from both ends inwards, while the sequencing direction of a long library is from the inside outwards. We can determine whether the directions of the paired-end sequences and their mapped congruent sequences are consistent based on the magnitude of the mapped coordinates of the start points of the sequences at both ends and the aforementioned direction (for a short library, the start point coordinate of the left end sequence is smaller than the end point coordinate, and the start point coordinate of the right end sequence is larger than the end point coordinate, indicating a consistent direction; for a long library, the opposite is true). Based on the direction determination results of the paired-end sequencing sequence and its mapped congruent sequence, determine whether the directions of two different congruent sequences are consistent. Calculate the distance between the start points of the two congruent sequences based on the direction determination results of the two different congruent sequences and the library length. Simultaneously, filter out ligations where one end is mapped to a position on the congruent sequence that exceeds the preset insertion length range, obtaining the ligation information from the sequencing sequence to each pair of congruent sequences. Figure 5 This is a schematic diagram illustrating the extraction of pairwise directional relationships and start-point distance information between consistent sequences from double-terminal sequence mapping data, as shown below. Figure 5 As shown, it specifically includes:

[0121] Define the stacking direction of the left end sequence mapping of the paired-end sequencing sequence as positive and its position as first. Assume the starting points of the stacking array on the left and right end mappings are F1 and F2, respectively, then:

[0122] a) If s1 < t1 and s2 < t2, then the directions of the two stacked array sequences are consistent with the paired-end sequencing sequences, that is, the directions of the two stacked arrays are opposite. Then the distance is: F2-F1 = in + s1 + s2, and the links len1-t1 > in + 3σ or len2-t2 > in + 3σ are screened out.

[0123] b) If s1 < t1 and s2 > t2, then the directions of the two stacked array sequences are opposite to those of the paired-end sequencing sequences, that is, the directions of the two stacked arrays are consistent. Then the distance is: F2 - F1 = in + s1 - s2, and the links len1 - t1 > in + 3σ or s2 > in + 3σ are eliminated.

[0124] c) If s1 > t1 and s2 > t2, then the directions of the two stacked array sequences are consistent with the paired-end sequencing sequences, that is, the directions of the two stacked arrays are opposite. Then the distance is: F2 - F1 = in - s1 - s2, and the links s1 > in + 3σ or s2 > in + 3σ are screened out.

[0125] d) If s1 > t1 and s2 < t2, then the directions of the two stacked array sequences are opposite to those of the paired-end sequencing sequences, that is, the directions of the two stacked arrays are consistent. The distance is: F2 - F1 = in - s1 + s2. Links that are s1 > in + 3σ or len2 - t2 > in + 3σ are eliminated.

[0126] If the input sequencing sequence is a single-molecule long-read sequencing sequence, then consider sequencing sequences mapped to multiple congruent sequences. First, assume that the coordinates of the alignment intervals aligned to multiple congruent sequences on the sequencing sequence, ordered by their starting points, are [qs1, qe1], [qs2, qe2], ..., [qs...]. k ,qe k The coordinates on the corresponding stacked matrix consistent sequence are [cs1, ce1], [cs2, ce2], ..., [cs k ,ce k Based on the above coordinates, the directional relationship and starting distance between any two congruent sequences aligned to the same sequencing sequence are used as the connection information for each pair of congruent sequences. Figure 6 This is a schematic diagram illustrating the directional relationships and origin distances between pairs of congruent sequences extracted from single-molecule long-read sequencing sequence mapping data.

[0127] Where the stacking array is aligned with the sequencing sequence orientation, then: F i =qs i -cs i If the stacking array is in the opposite direction to the sequencing sequence, then: F i =qs i -cs i +len i The distance is: F j -F i , of which F i and F j The starting point of the overlapping sequences on the sequencing sequence mapping; finally, the overlapping links on the same sequencing sequence are screened out.

[0128] S14: Construct an undirected weighted overlay matrix based on the uniquely mapped sequencing sequence information and the selected consistent sequence set, and divide the undirected weighted overlay matrix into multiple connected components, including:

[0129] S1401: Defines the threshold parameter m for the number of connections between stacked arrays. The setting of the threshold parameter mainly depends on the depth of the sequencing data. It can be designed to be half of the total depth of the sequencing data, and should not be lower than 3.

[0130] S1402: Using the consistent sequence set of the assembled stacked array as vertices, an edge is formed between the single-mapped sequencing sequence stacked arrays if the number of connections between them is greater than the threshold parameter m, thus constructing an undirected weighted stacked array graph; this edge corresponding to n connections contains the direction and distance information between the two vertices, with some connections indicating that the two vertices have the same direction and the other part indicating that their directions are opposite.

[0131] S1403: Let a denote the number of same-direction connections and the number of non-direction connections between vertices i and j in an undirected weighted matrix. ij and b ij The corresponding distance information is entered into the stacking distance lists of the same direction and opposite direction, respectively;

[0132] S1404: Divide the undirected weighted overlay matrix into multiple connected components.

[0133] S15: The direction of the stacking matrix corresponding to each vertex in each connected component is determined using a combinatorial optimization model, including:

[0134] S1501: with a ij and b ij The maximum value in the matrix assigns a weight to each edge in the connected component and marks the direction of that edge. If a ij ≥b ij Then the weight w ij =a ij D ij =1; otherwise, weight w ij =b ij D ij =-1; Traverse all vertices to obtain the initial direction of each vertex in the connected component;

[0135] S1502: Iteratively optimize the total number of conflicting directional connections on the current connected component. In each iteration, calculate the single-point change of the objective function caused by changing each vertex. If there is a vertex whose single-point change is negative, select the vertex with the largest ratio of single-point change to total number of connections and change its direction. After updating the direction, continue to the next iteration until changing any vertex can no longer reduce the total number of conflicts. Stop the iteration and obtain the direction of the stack matrix corresponding to each vertex. Figure 7 The algorithm flow for optimizing vertex orientation assignments using a combinatorial optimization model on a connected component of a stacked matrix is ​​shown.

[0136] In some embodiments, a depth-first traversal method with decreasing weights is used to traverse all vertices to obtain the initial direction of each vertex in the connected component: starting from vertex zero (direction set to 1), always visiting the current vertex (the meaning of visiting: visiting vertex j from vertex i, let D...). j =D ij D jThen, select the neighbor with the largest weight of the current vertex and iterate through all vertices to obtain the initial direction of all vertices.

[0137] In some other embodiments, the maximum spanning tree method is used to traverse all vertices to obtain the initial direction of each vertex of the connected component: first, find a maximum weight spanning tree for the current connected component, and then traverse the connected component using the edges of the spanning tree to obtain the initial direction of all vertices.

[0138] S16: Establish a regression model to estimate the arrangement coordinates of the stacked matrices on the genome being tested. Sort all stacked matrices in the same connected component according to the orientation of the stacked matrices corresponding to each vertex and the arrangement coordinates to obtain the quasi-architecture of that connected component. Obtain an architecture set based on the quasi-architectures of multiple connected components, including:

[0139] S1601: Establish a regression model;

[0140] S1602: Based on the connection information of the pairwise overlapping matrix sequence and the direction of the overlapping matrix corresponding to each vertex, determine the type of each edge, including edges in the same direction and edges in opposite directions;

[0141] S1603: Based on the type of each edge, extract the median of the stacked matrix distance list corresponding to each edge as the observation value of the regression model;

[0142] S1604: Set the maximum error threshold for the residuals of the weighted truncated least squares estimation. To ensure the accuracy of robust regression, we set this threshold as small as possible, but considering the variance of the library length and the efficiency of iteration, we set it to 100.

[0143] S1605: The initial weight is the product of the inverse of the variance estimate of the length of the pool and the number of connections in the current edge that conform to the direction of the current vertex. Calculate the weighted least squares estimate of the starting coordinates of the consistent sequence of stacked matrices. In the coefficient matrix X, the first row is (1,0,...,0), and the remaining rows represent the edges between the corresponding vertices i and j. Therefore, the i-th item is -1 and the j-th item is 1.

[0144] S1606: Calculate the residuals of all samples based on the previous estimate. Stop the iteration when the maximum residual value is less than the residual allowable threshold `maxError`, and obtain a feasible estimate. Otherwise, select samples whose residuals are less than the (100-iter) percentile, where iter represents the current iteration number;

[0145] Determine if the coefficient matrix corresponding to the selected samples is of full rank: if not, further divide the current subgraph into several connected components based on the coefficient matrix, and restart from step S1501 on each connected component of the subgraph; if it is of full rank, perform weighted least squares estimation again on the selected sample subset to obtain... Repeat step S1606 until a feasible estimate is obtained.

[0146] S1607: Sort all the stacked matrices in the current connected component according to the direction of the stacked matrices corresponding to each vertex, and organize them into the same coordinate system to obtain a quasi-architecture. Obtain an architecture set based on the quasi-architectures of multiple connected components.

[0147] S17: The coordinates of the stacked arrays for each architecture in the architecture set are sorted out to eliminate overlapping or redundant arrays, resulting in the assembled genome architecture, including:

[0148] S1701: Adjust all the consistent sequences of stacked matrices in the architecture according to the direction of the stacked matrices corresponding to each vertex, swap the starting and ending coordinates of the reversed stacked matrices, so that the starting coordinates of all stacked matrices in the architecture arrangement are before the ending coordinates, and sort the stacked matrices in ascending order of starting coordinates.

[0149] S1702: Starting from the first stacked matrix after sorting, proceed sequentially, performing local sequence alignment on the stacked matrices that are contained in or intersected by the last selected stacked matrix at the current coordinates. If the alignment is successful, remove the contained stacked matrix, or merge the intersected stacked matrix into the current stacked matrix; select stacked matrices that are adjacent to the endpoint coordinates of the current stacked matrix and have connection information to add them to the current architecture, and leave stacked matrices that are adjacent to the endpoint coordinates of the current stacked matrix but have no connection information; where...

[0150] 1) If the coordinates of two overlapping matrices i and j satisfy:

[0151] 0 < start j -end i ≤2kbp, and end j >end i Then the two stacked arrays are adjacent;

[0152] 2) If the coordinates of two identical sequences i and j of two superimposed matrices satisfy:

[0153] start i <start j <end j <end i Then i contains j;

[0154] 3) If the coordinates of two identical stacked matrices, i and j, satisfy:

[0155] start i <start j <end i <end j Then the two overlapping matrices intersect.

[0156] Where start and end represent the starting and ending coordinates of the stacked matrix consistent sequence in the quasi-architecture, respectively;

[0157] S1703: Update the current stack to the latest stack added to the current architecture, repeat step S1702 until the last stack on the architecture is reached, end the current architecture, and output to the final assembled architecture file;

[0158] S1704: Repeat steps S1702-S1703 in the retained stacks to determine the next architecture arrangement and output it to the final assembled architecture file; repeat this step until all stacks are arranged into a final architecture, to obtain the assembled genome architecture, the position and orientation information of the stacks in the architecture, the estimated length of the gap between adjacent stacks and the estimated variance, and the new stack consistent sequence obtained by merging the coordinate cross and successful alignment.

[0159] Figure 9 This diagram illustrates the final architectural arrangement, utilizing regression-estimated coordinates, connectivity information, and local comparisons of overlapping regions.

[0160] In some other embodiments, to obtain a longer architecture, the quasi-architecture set comprising the quasi-architectures is replaced by the set of consistent sequences of overlapping matrices, and steps S14-S16 are repeated to obtain the architecture set. During the iteration process, the following should be noted:

[0161] A. The connections in the stacking diagram come directly from the information of the sequencing sequence aligned to two different stacking diagrams with the same sequence. In the quasi-architecture diagram, it is necessary to select the sequencing sequence information aligned to two different quasi-architectures, that is, to select the mapping information of two stacking diagrams with the same sequence that are not in the same quasi-architecture.

[0162] B. Updated from the definition of connection in A, the connection between two quasi-archites, indicating whether they are in the same or opposite direction, needs to supplement the consideration of the direction of the consistent sequences of the overlay matrix mapped at both ends in their respective quasi-archites. If the consistent sequences of the overlay matrix mapped at both ends are in the same direction and have the same direction in their respective quasi-archites, then the directions of the two quasi-archites are the same; if the consistent sequences of the overlay matrix mapped at both ends are in opposite directions and have opposite directions in their respective quasi-archites, then the directions of the two quasi-archites are also the same; otherwise, the directions of the two quasi-archites are opposite. Due to the different connection definition in A, the formula for calculating the distance between the starting points of the two quasi-architecture sequences also needs to be replaced according to different cases, which are divided into the following 8 cases, where start and end represent the coordinates of the starting point and ending point of the consistent sequence of the overlay matrix in the quasi-architecture, respectively, and diff is the distance between the starting points of the consistent sequence of the overlay matrix calculated in step S1303:

[0163] a) If the overlapping matrices mapped at both ends are of opposite direction and have opposite directions in their respective quasi-architectures, with the former being positive and the latter negative, then:

[0164]

[0165] b) If the overlapping matrices mapped at both ends have opposite directions but the same positive direction in their respective quasi-architectures, then:

[0166]

[0167] c) If the consistent sequences of the overlapping matrices mapped at both ends are in the same direction and have the same positive direction in their respective quasi-architectures, then:

[0168]

[0169] d) If the consistent sequences of the overlapping matrices mapped at both ends are in the same direction but have opposite directions in their respective quasi-architectures, with the former being positive and the latter negative, then:

[0170]

[0171] e) If the consistent sequences of the overlapping matrices mapped at both ends are in the same direction and have the same direction in their respective quasi-architectures, both being negative, then:

[0172]

[0173] f) If the consistent sequences of the overlapping matrices mapped at both ends are in the same direction but have opposite directions in their respective quasi-architectures, with the former being negative and the latter positive, then:

[0174]

[0175] g) If the overlapping matrices mapped at both ends are of opposite direction and have opposite directions in their respective quasi-architectures, with the former being negative and the latter positive, then:

[0176]

[0177] h) If the overlapping matrices mapped at both ends have identical sequences that are out of direction but have the same direction in their respective quasi-architectures, both being negative, then:

[0178]

[0179] Figures 8a-8b This is a schematic diagram illustrating eight scenarios for calculating the pairwise directional relationships and distance information between quasi-archites, such as... Figures 8a-8b The diagram shows the orientation of the overlapping sequence mapped to both ends of the paired-end sequencing sequence in their respective quasi-architectures.

[0180] C. Appropriate changes should also be made in the selection of parameters. The threshold parameter m' for the number of connections between vertices of the quasi-architecture should be 1.5 times m, and the allowable threshold maxError′ for the residuals of the weighted truncated least squares estimation should be 1.5 times maxError.

[0181] A fourth aspect of the present invention provides a genome architecture assembly system, which performs genome architecture assembly using the described genome architecture assembly method.

[0182] The genome architecture assembly method provided by this invention employs a hierarchical iterative algorithm to fully utilize all valuable sequencing information and adopts a strict-then-relaxed strategy to ensure the accuracy of the underlying quasi-architecture before iterative extension. The regression-based array arrangement calculation can also provide the error of the gap length estimation between adjacent arrays, that is, it provides an interval estimate, providing a reliable reference interval for subsequent gap filling. It can simultaneously utilize second-generation paired-end sequencing and third-generation single-molecule long-read sequencing, or it can use only one type of sequencing data, showing strong adaptability to the platform.

[0183] This invention also provides a method for assembling genome sequences, such as... Figure 10 As shown, the method includes:

[0184] S21: Extract a third sequencing sequence set at a certain depth, which is used to generate an assembled stacked array set according to the aforementioned genome stacked array assembly method. The assembled stacked array set is used to generate an assembled genome architecture according to the aforementioned genome architecture assembly method.

[0185] S22: Extract a fourth sequencing sequence set from another depth to update and improve the assembled genome architecture, and perform a quality assessment on the improved assembled genome. If the quality meets the standard, the improved assembled genome is stored as a candidate assembled genome, including:

[0186] S2201: The fourth sequencing sequence The sequence is mapped to the assembled genome, and sequencing sequences that are mapped to multiple locations are removed to form a stack of mono-mapped sequencing sequences. For each column in the stack, if it is covered by more than d0 mono-mapped sequences, the most frequent base in that column is taken to update the corresponding base in the assembled genome, resulting in the improved assembled genome.

[0187] S2202: The fourth sequencing sequence The genome is mapped to the improved assembled genome, and its quality is evaluated based on the mapping results. If the quality meets the standard, it is stored as a candidate assembled genome.

[0188] S23: Repeat steps S21 and S22 until a predetermined number of candidate assembled genomes are obtained. Perform multiple sequence alignment on all candidate assembled genomes and connect the most frequent bases in each aligned column to obtain the final assembled genome.

[0189] The final assembled genome is stored as a file in a standard format.

[0190] The present invention also provides a genome sequence assembly system, which uses the aforementioned genome sequence assembly method to assemble genome sequences.

[0191] This invention achieves genome stack assembly based on the regression concept, and also achieves genome architecture assembly based on the regression concept. Finally, the sequencing data is independently sampled, assembled, improved, and evaluated multiple times through resampling technology, and then integrated into the final assembled genome. This can reduce the uncertainty of the assembly results caused by noise in the sequencing data, and reduce the sensitivity of the assembly results to the selection of assembly parameters.

[0192] This invention also provides a machine-readable storage medium storing computer program instructions that, when executed by a processor, implement the above-described genome stack assembly method.

[0193] Embodiments of the present invention also provide another machine-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the above-described genome architecture assembly method.

[0194] This invention also provides a machine-readable storage medium storing computer program instructions that, when executed by a processor, implement the above-described genome sequence assembly method.

[0195] Those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. This program is stored in a storage medium and includes several instructions to cause a microcontroller, chip, or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as a USB flash drive, a portable hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0196] The optional embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the embodiments of the present invention are not limited to the specific details described above. Within the scope of the technical concept of the embodiments of the present invention, various simple modifications can be made to the technical solutions of the embodiments of the present invention, and these simple modifications all fall within the protection scope of the embodiments of the present invention. It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the embodiments of the present invention will not further describe the various possible combinations.

[0197] Furthermore, various different embodiments of the present invention can be combined in any way, as long as they do not violate the spirit of the embodiments of the present invention, they should also be regarded as the content disclosed by the embodiments of the present invention.

Claims

1. A method for assembling a genome stack, characterized in that, The method includes: S01: Extract the first sequencing sequence set at a certain depth; S02: The sequencing sequences in the first sequencing sequence set are compared pairwise using a sequence alignment algorithm to obtain pairwise alignment information of the sequencing sequences; S03: Based on the pairwise alignment information of the sequencing sequences, all alignments that meet the overlap condition are collectively represented as the first linear regression model; S04: Solve the first linear regression model using an iterative reweighted least squares algorithm to obtain the solution with the minimum global loss function and a robust estimate of the coordinates of the sequencing sequences in the first sequencing sequence set on the genome being tested, including: S041: Solving sparse linear equations about β using an accelerated algorithm for large sparse linear equation systems: The solution obtained is as The initial value; S042: In each iteration, based on the current... Calculate the residual vector ,as well as Maximum value of the absolute value of each component ; S043: If Then calculate the weight matrix. ; in, The preset convergence threshold is used; the weighting function is... ; S044: Based on the calculated W, assign the smaller of the weight values ​​of the two samples that are each other's dual to these two samples; S045: Solving sparse linear equations about β using an accelerated algorithm for large sparse linear equation systems: ; S046: Update using the obtained solution Delete the rows in matrices X and Y with weights of 0, and return to the next iteration; S047: Until Stop the iteration and record the matrix X of the last iteration; S048: For the final round Output after rounding each component This is a robust estimate of the coordinates of the sequencing sequences in the first sequencing sequence set on the genome being tested; S05: The sequencing sequence is divided into multiple connected components by an undirected graph constructed based on the solution that minimizes the global loss function. The sequencing sequences are then arranged according to a robust estimate of the coordinates of the sequencing sequences within each connected component to generate a stacked matrix. S06: Merge the stacked arrays based on the pairwise alignment information of the sequencing sequences to obtain a preliminary assembled stacked array set; S07: Obtain the consistent sequence of each array in the preliminary assembled array set, map the sequencing sequence to the consistent sequence, and obtain the mapping result of the sequencing sequence to the consistent sequence. The mapping result is the new assembled array set.

2. The genome array assembly method according to claim 1, characterized in that, Step S03 includes: S031: Assuming the sequencing sequences in the first sequencing sequence set... Align the rightmost base to the sequencing sequence internal, The right overhang contains One base, thus obtaining an observation sample: ; in, yes The coordinates on the single strand of DNA it is located at. yes The coordinates on the single strand of DNA it is located at. Indicates the observation error of the observed sample; S032: Assuming a second sequencing sequence Align the leftmost base to the first sequencing sequence internal, The left overhang contains Each base pair yields a dual observation sample: ; in, yes The coordinates of the inverse complementary sequence on the complementary single chain, yes The coordinates of the inverse complementary sequence on the complementary single chain, This represents the observation error of the dual observation sample; S033: Based on the pairwise alignment information of the sequencing sequences, all aligned observation samples and their pairs that meet the overlap condition are integrated into a matrix equation: ; Each row represents an observation sample, and the vector... For the parameters to be estimated, X is a large sparse matrix. Each row of X has only two non-zero elements, 1 and -1. X is stored in sparse matrix format.

3. The genome array assembly method according to claim 1, characterized in that, Step S05 includes: S051: Using the sequencing sequences in the first sequencing sequence set and the reverse complementary sequences corresponding to the sequencing sequences as vertices, determine two sequencing sequences based on the two non-zero elements in each row of the matrix X in the last round of iteration, establish an edge between the vertices corresponding to these two sequencing sequences, and construct an undirected graph. S052: Traverse the undirected graph and divide it into multiple connected components, obtaining 2k connected components, denoted as . ,in and The sequencing sequences contained therein correspond one-to-one and are reverse complementary sequences to each other; S053: Arrange the sequencing sequences in each connected component according to a robust estimate of their coordinates on the genome being sequenced to obtain the stacked matrix set: ; Stacked array and These correspond to two subsequences on the DNA double helix that are inversely complementary.

4. The genome array assembly method according to claim 3, characterized in that, Step S06 includes: S061: Use the stacked array in the stacked array set as the vertex; S062: For any two superimposed matrices, if The 3' end s sequencing sequences and If s sequencing sequences at the 5' end have more than t alignments that satisfy the overlap condition, then a sequence is established. Corresponding vertex pointing Construct a directed graph using the directed edges corresponding to the vertices; where the weight of the edge is set to the number of logs that satisfy the stacking condition, and the values ​​of s and t are positively correlated with the sampling depth; S063: The stacked array containing the most sequencing sequences from the directed graph. Begin by determining the stacked array The stacked matrix pointed to by the edge with the largest weight at the starting point ; S064: From Select a sequencing sequence from the 3' end , and from Select a sequencing sequence from the 5' end. , and The comparison between them satisfies the overlap condition; S065: Assumption Align the rightmost base to internal, The right overhang contains If there are 1 base, then... Robust estimates of the coordinates of all sequencing sequences are added ; S067: Will and The sequencing sequences in the array are arranged and merged into a stacked array according to the robust estimates of their coordinates, while simultaneously... and Remove the corresponding vertex in the directed graph; S068: From Repeat steps S063-S067 until no more stacks can be merged. Then, starting from the stack with the most sequencing sequences among the remaining vertices in the directed graph, repeat steps S063-S067 until all vertices in the directed graph have been traversed. The resulting merged stacks are taken as the initial assembled stack set, denoted as . .

5. The genome array assembly method according to claim 4, characterized in that, Step S07 includes: S071: For any one of the stacked arrays in the preliminary assembled stacked array set Take the base with the highest frequency in each column to form a stacked array. The corresponding consistent sequence is denoted as ; S072: For stacked arrays Any sequencing sequence The leftmost base is at The coordinates in the middle are as follows , The rightmost base is at The coordinates in the middle are as follows ,Will The Middle From the first base to the second Extracting the subsequence between bases and... Perform local sequence alignment; S073: The portion of the local sequence that matches the alignment is taken as... Towards The mapping results All sequencing sequences to The mapping results constitute a new superposition matrix. ; For the updated assembled array set.

6. A method for assembling a genome architecture, characterized in that, The method includes: S11: Extract a second sequencing sequence set at a certain depth; S12: Assemble the genome stack array according to the second sequencing sequence set using the genome stack array assembly method of any one of claims 1-5 to obtain the assembled stack array set. ; S13: Combine the second sequencing sequence set with the assembled stack set The sequence is compared with the consistent sequence set. Based on the comparison results, the single-mapped sequencing sequences and their mappings are selected. The single-mapped sequencing sequences and their mappings are then selected according to the set conditions to obtain the uniquely mapped sequencing sequence information and the selected consistent sequence set. S14: Construct an undirected weighted overlapping matrix based on the uniquely mapped sequencing sequence information and the selected consistent sequence set, and divide the undirected weighted overlapping matrix into multiple connected components; S15: Use a combinatorial optimization model to determine the direction of the stacking matrix corresponding to each vertex in each connected component; S16: Establish a regression model to estimate the arrangement coordinates of the stacked arrays on the genome being tested. Sort all the stacked arrays in the same connected component according to the direction of the stacked arrays corresponding to each vertex and the arrangement coordinates to obtain the quasi-architecture of the connected component. Obtain the architecture set based on the quasi-architectures of multiple connected components. S17: The coordinates of the stacked array of each architecture in the architecture set are sorted out to eliminate the contradiction of coordinate overlap or redundant stacked arrays, so as to obtain the assembled genome architecture.

7. The genome architecture assembly method according to claim 6, characterized in that, The step of selecting single-mapped sequencing sequence stacks based on alignment results includes: S1301: For paired-end sequencing sequences, the mapping results from the alignment results to the two ends of the paired-end sequence onto the same congruent sequence are selected to estimate the library length of the paired-end sequencing fragment library. The mapping results include the coordinates of the starting point of each end sequence alignment onto the congruent sequence, denoted as s1,t1 and s2,t2 respectively; the library length observation is denoted as: For each fragment library, a predetermined number of observations are taken, and the median M of these observations is used as the observation for the library length. The standard deviation of this median is then calculated. ; S1302: The mapping coordinates of the paired-end sequencing sequences mapped to two different stacked consistent sequences are still denoted as s1,t1 and s2,t2, respectively. S1303: Determine whether the orientation of the paired-end sequencing sequence and its mapped contiguous sequence is consistent based on the mapping coordinates; determine whether the orientations of two different contiguous sequences are consistent based on the orientation determination results of the paired-end sequencing sequence and its mapped contiguous sequence; calculate the distance between the starting points of the two contiguous sequences based on the orientation determination results of the two different contiguous sequences and the library length; filter out the connections where one end is mapped to the contiguous sequence and the position exceeds the preset insertion length range, and obtain the connection information from the sequencing sequence to the pairwise contiguous sequences.

8. The genome architecture assembly method according to claim 7, characterized in that, Step S1303 includes: The stacking direction on the left end sequence mapping of the paired-end sequencing sequence is defined as positive and its position is first. It is assumed that the starting points of the stacking arrays on the left and right end mappings are respectively... , ,but: a) If and If the orientations of the two congruent sequences in the two stacked arrays are consistent with those of the paired-end sequencing sequences, meaning the orientations of the two stacked arrays are opposite, then the distance is: , screening or The connection; b) If and If the orientations of the two congruent sequences in the two stacked arrays are opposite to those of the paired-end sequencing sequences, then the distance between them is: , screening or The connection; c) If and If the orientations of the two congruent sequences in the two stacked arrays are consistent with those of the paired-end sequencing sequences, meaning the orientations of the two stacked arrays are opposite, then the distance is: , screening or The connection; d) If and If the orientations of the two congruent sequences in the two stacked arrays are opposite to those of the paired-end sequencing sequences, then the distance between them is: , screening or The connection.

9. The genome architecture assembly method according to claim 6, characterized in that, The step of selecting single-mapped sequencing sequences and their mappings based on alignment results includes: S1311: For single-molecule long-read sequencing sequences, it is assumed that the coordinates of the alignment intervals aligned to multiple congruent sequences on the sequencing sequence, ordered by the starting point, are as follows: The coordinates on the corresponding stacked matrix consistent sequence are as follows: Based on the above coordinates, the directional relationship and starting distance between any two congruent sequences aligned to the same sequencing sequence are used as the connection information for each pair of congruent sequences. If the stacking array is aligned with the sequencing sequence orientation, then: If the stacking array is in the opposite direction to the sequencing sequence, then: The distance is: ; S1312: Screening out overlapping links on the same sequencing sequence.

10. The genome architecture assembly method according to claim 7 or 9, characterized in that, Step S14 includes: S1401: Defines the threshold parameter m for the number of connections between stacked matrices to form edges; S1402: Using the consistent sequence set of the assembled stacked array as the vertex, an edge is formed when the number of connections between the stacked arrays of single-mapped sequencing sequences is greater than the threshold parameter m, and an undirected weighted stacked array graph is constructed. S1403: Let the number of same-direction connections and the number of opposite-direction connections between vertices i and j in an undirected weighted overlapping matrix be denoted as follows: and The corresponding distance information is entered into the stacking distance lists of the same direction and opposite direction, respectively; S1404: Divide the undirected weighted overlay matrix into multiple connected components.

11. The genome architecture assembly method according to claim 10, characterized in that, Step S15 includes: S1501: with and The maximum value in the matrix assigns a weight to each edge in the connected component and marks the direction of that edge. Then the weight , Otherwise, the weight , ; Traverse all vertices to obtain the initial direction of each vertex in the connected component; S1502: Iteratively optimize the total number of conflicting directional connections on the current connected component. In each iteration, calculate the single-point change of the objective function caused by changing each vertex. If there is a vertex with a negative single-point change, select the vertex with the largest ratio of single-point change to total number of connections and change its direction. After updating the direction, continue to the next iteration until changing any vertex can no longer reduce the total number of conflicts. Stop the iteration and obtain the direction of the stack matrix corresponding to each vertex.

12. The genome architecture assembly method according to claim 11, characterized in that, Step S16 includes: S1601: Establish a regression model; S1602: Determine the type of each edge based on the connection information of the pairwise overlapping matrix sequence and the direction of the overlapping matrix corresponding to each vertex; S1603: Based on the type of each edge, extract the median of the stacked matrix distance list corresponding to each edge as the observation value of the regression model; S1604: Set the allowable threshold for the residuals of weighted truncated least squares estimates. ; S1605: The initial weight is the product of the inverse of the variance estimate of the length of the pool and the number of connections in the current edge that conform to the direction of the current vertex. Calculate the weighted least squares estimate of the starting coordinates of the consistent sequence of stacked matrices. The coefficient matrix is ​​used for this purpose. Except for the first row (1,0,...,0), the remaining rows represent the edges between the corresponding vertices i and j. Therefore, the i-th item is -1 and the j-th item is 1. S1606: Calculate the residuals of all samples based on the previous round of estimation. When the maximum residual value is less than the residual allowable threshold... Stop iterating to obtain a feasible estimate. Otherwise, select samples with residuals less than A sample of quantiles, where, Indicates the current iteration number; Determine if the coefficient matrix corresponding to the selected samples is of full rank: if not, further divide the current subgraph into several connected components based on the coefficient matrix, and restart from step S1501 on each connected component of the subgraph; if it is of full rank, perform weighted least squares estimation again on the selected sample subset to obtain... Repeat step S1606 until a feasible estimate is obtained. ; S1607: Sort all the stacked matrices in the current connected component according to their corresponding feasible estimates based on the direction of the stacked matrices corresponding to each vertex, organize them into the same coordinate system, and obtain a quasi-architecture. Obtain an architecture set based on the quasi-architectures of multiple connected components.

13. The genome architecture assembly method according to claim 12, characterized in that, Step S17 includes: S1701: Adjust all the consistent sequences of stacked matrices in the architecture according to the direction of the stacked matrices corresponding to each vertex, swap the starting and ending coordinates of the reversed stacked matrices, so that the starting coordinates of all stacked matrices in the architecture arrangement are before the ending coordinates, and sort the stacked matrices in ascending order of starting coordinates. S1702: Starting from the first stacked matrix after sorting, proceed sequentially to the next stacked matrix. Perform local sequence alignment on the stacked matrices that are contained in or intersected by the last stacked matrix currently selected. If the alignment is successful, remove the contained stacked matrix or merge the intersected stacked matrix into the current stacked matrix. Select the stacked matrix that is adjacent to the current stacked matrix's endpoint coordinates and has connection information, and add it to the current architecture; leave the stacked matrix that is adjacent to the current stacked matrix's endpoint coordinates but has no connection information; among them... 1) If the coordinates of two overlapping matrices i and j satisfy: ,and Then the two stacked arrays are adjacent; 2) If the coordinates of two overlapping matrices i and j satisfy: Then i contains j; 3) If the coordinates of two overlapping matrices i and j satisfy: Then the two overlapping matrices intersect. in, These represent the starting and ending coordinates of the consistent sequence in the quasi-architecture, respectively. S1703: Update the current stack to the latest stack added to the current architecture, repeat step S1702 until the last stack on the architecture is reached, end the current architecture, and output to the final assembled architecture file; S1704: Repeat steps S1702-S1703 in the retained stacks to determine the next architecture arrangement and output it to the final assembled architecture file; repeat this step until all stacks are arranged into a final architecture, to obtain the assembled genome architecture, the position and orientation information of the stacks in the architecture, the estimated length of the gap between adjacent stacks and the estimated variance, and the new stack consistent sequence obtained by merging the coordinate cross and successful alignment.

14. The genome architecture assembly method according to claim 12, characterized in that, The method further includes: replacing the set of consistent sequences of stacked matrices with the set of quasi-architectures composed of the quasi-architectures, and repeating steps S14-S16 to obtain the set of architectures.

15. A method for assembling a genome sequence, characterized in that, The method includes: S21: Extract a third sequencing sequence set at a certain depth for use in generating an assembled stacked array set using the genome stacked array assembly method according to any one of claims 1-5, wherein the assembled stacked array set is used in generating an assembled genome architecture using the genome architecture assembly method according to any one of claims 6-14. S22: Extract a fourth sequencing sequence set from another depth to update and improve the assembled genome architecture, and evaluate the quality of the improved assembled genome. If the quality meets the standard, the improved assembled genome is stored as a candidate assembled genome. S23: Repeat steps S21-S22 until a predetermined number of candidate assembled genomes are obtained. Perform multiple sequence alignment on all candidate assembled genomes and connect the most frequent bases in each aligned column to obtain the final assembled genome.

16. The genome sequence assembly method according to claim 15, characterized in that, Step S22 includes: S2201: The fourth sequencing sequence Mapping to the assembled genome, removing sequencing sequences that simultaneously map to multiple locations, forming a mono-mapped sequencing sequence stack; for each column in the stack, if it is mapped to more than... If a single-image sequence covers the entire sequence, the most frequent base in that sequence is used to update the corresponding base position in the assembled genome, resulting in an improved assembled genome. S2202: The fourth sequencing sequence The genome is mapped to the improved assembled genome, and its quality is evaluated based on the mapping results. If the quality meets the standard, it is stored as a candidate assembled genome.

17. A genome stacking assembly system, characterized in that, The system uses the genome stacking assembly method according to any one of claims 1-5 to perform genome stacking assembly.

18. A genome architecture assembly system, characterized in that, The system uses the genome architecture assembly method according to any one of claims 6-14 to assemble the genome architecture.

19. A genome sequence assembly system, characterized in that, The system uses the genome sequence assembly method described in claim 15 or 16 to assemble genome sequences.

Images