A second-generation de novo assembly method and system based on gene numerical expression

Abstract

The present invention discloses a second-generation de novo assembly method and system based on digital gene expression. The method comprises the following steps: S1: implementing base sequencing data management through a flexible disk array (RAID); S2: gene base sequence analysis and custom numbering; S3: calculating the front and back centroid values ​​of two groups of specific lengths for each sequencing read; S4: constructing a B+ tree index based on the base sequence and read ID of each sequencing read; S5: generating an ID docking table using an artificial intelligence algorithm or a numerical fast sort algorithm; S6: multi-threaded assembly of small contigs based on the ID docking table, outputting a thread identification table, small contig fragments and their numbers, and a single nucleotide variant (SNP) information table; S7: constructing a De Brujin Graph based on the thread identification and pathfinding to assemble the small contigs into contigs; and S8: assembling a scaffold sequence based on the contigs. The present invention completely solves the problem of excessive memory usage in bioinformatics assembly and reduces the risk of assembly errors.

Description

Technical Field

[0001] The present invention belongs to the technical field of bioinformatics, and in particular relates to a second-generation de novo assembly method and system based on numerical gene expression. Background Art

[0002] Since Roche launched its first second-generation sequencer, the Roche 454, in 2005, life sciences have officially entered the era of high-throughput sequencing. The introduction of Illumina's sequencing platforms has significantly reduced the price of second-generation sequencing, making high-throughput sequencing widely available across various life science research fields. To this day, second-generation short-read sequencing technology maintains a dominant position in the global sequencing market.

[0003] In the field of bioinformatics, genome assembly is a core issue in genome analysis. Although second-generation sequencing technology has greatly increased sequencing throughput, the individual sequences it obtains are very short, often only 50 to 300 base pairs. Furthermore, due to the deeper sequencing coverage, the number of sequences that need to be processed in genome projects has increased rapidly. These shorter individual sequences need to be spliced ​​into longer, complete sequences for further analysis to improve the accuracy of species annotation analysis.

[0004] Due to the limitations of current assembly technologies and the complexity of real-world scenarios, the final assembled sequence may still differ from the true genome sequence, or even result in only a few sequences that cannot be further connected. The current mainstream assembly algorithm is the De Brujin Graph (DBG) assembly algorithm. Since reads (sequencing fragments) generated by second-generation sequencing are short and numerous, it is difficult to detect overlap between reads. Therefore, reads are converted into fixed-length k-mers (nucleotide sequences of length K) to identify overlaps and construct a De Bruijn graph, transforming the short sequence assembly problem into an Euler path problem within the De Bruijn graph. However, the k-mer-based DBG assembly algorithm still has significant drawbacks, such as the short kmer sequences are susceptible to sequencing errors, and the inability to perform fast sorting when finding overlaps between bases. Summary of the Invention

[0005] The purpose of the present invention is to provide a second-generation de novo assembly method and system based on digital gene expression to solve the technical problems of slow sequencing read arrangement speed, excessive memory of bioinformatics assembly, and high risk of assembly errors in genome assembly.

[0006] The purpose of the present invention is to achieve the following technical solution: a second-generation de novo assembly method based on gene digital expression, comprising the following steps:

[0007] S1: Base sequencing data management is achieved through soft disk array RAID;

[0008] S2: Gene analysis base sequence and custom numbering;

[0009] S3: Calculate the front and back centroid values ​​of two groups of specific lengths for each sequencing fragment read;

[0010] S4: Construct a B+ tree index based on the base sequence and read ID of each sequencing fragment read;

[0011] S5: Apply an artificial intelligence algorithm to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical rapid sorting algorithm to sort the centroid values ​​in sequence, and generate an ID docking table based on the rapid sorting results;

[0012] S6: Multi-threaded assembly of small contigs based on the ID docking table, outputting a thread identification table, small contig fragments and their numbers, and a single nucleotide variant (SNP) information table;

[0013] S7: Construct De Brujin Graph based on thread identifiers and find paths, assemble small contigs into contigs;

[0014] S8: Assemble scaffolds based on contigs.

[0015] Furthermore, step S1 specifically includes: through a built-in RAID card or NFS protocol, the base sequencing data stored in the disk is gradually transferred into a pseudo-virtual memory composed of a high IOPS hard disk in the form of a data stream, and stored in a dispersed manner, wherein the base sequencing data includes centroid data and gene files.

[0016] Furthermore, step S2 specifically includes: reading the gene file from the solid-state drive with high random read and write capabilities through the CPU to obtain multiple gene sequencing fragment reads, and numbering each gene sequencing fragment read to obtain a random number ID, and then transmitting the gene sequencing fragment reads and the random number ID to the GPU through the CUDA stream.

[0017] Furthermore, the step S3 specifically includes: using the GPU to perform macro-layer identification on the base sequencing data, removing m and n bases before and after each sequencing fragment read, respectively, where the values ​​of m and n are different, to obtain two groups of four split front and back layer fragments, and using unsupervised machine learning to express the two groups of front and back layer sequences as centroid values ​​calculated by hard clustering.

[0018] Furthermore, the step S4 specifically includes: the base sequence of each sequencing fragment read and its corresponding Should The numbered read IDs are used to construct a B+ tree index to achieve random extraction of base sequences.

[0019] Furthermore, step S5 specifically includes: after obtaining the front and back centroid values ​​of two groups of different fragments of all sequencing sequences, using an artificial intelligence algorithm, taking the two groups of front and back layer centroids and sequence lengths of the two sequences as input vectors, and whether the two sequences can match as the output vector, applying a deep learning neural network (DNN) to obtain a matching model of the relationship between the front and back centroid values ​​and the base sequence alignment, and generating a highly accurate ID docking table based on this, or using a numerical quick sort algorithm to sort the centroids in sequence, and generating an ID docking table based on this.

[0020] Furthermore, the artificial intelligence algorithm includes the following sub-steps:

[0021] S51: transform the two sets of front and back centroid values ​​of the two sequences and their sequence lengths to within [-1, 1] as the input vector x of the data set;

[0022] S52: Use the KMP algorithm to align the two sequences. If they can be aligned, it returns TRUE, y=1; if they cannot be aligned, it returns FALSE, y=0, and y is used as the output vector of the data set.

[0023] S53: Before training, first set a random distribution with a mean of 0 to initialize the network weights. During training, the input vector x is first sent to the input layer of the DNN network for forward propagation. After passing through several hidden layers, the forward output y' is obtained;

[0024] S54: Calculate the error between y' and the true value y, and then backpropagate the error from the output layer to the first hidden layer, adjusting the connection weights of the network in the direction of reducing the error;

[0025] S55: The above weight adjustment process is continuously repeated until the training error converges to a minimum value, that is, the DNN network weights that meet the requirements are obtained.

[0026] Furthermore, step S6 specifically includes: multi-threaded reading of the numbered read ID in the ID docking table, using a B+ tree index to extract the corresponding numbered read ID from the hard disk or virtual memory, and matching them with the benchmark sequencing fragment reads in turn, matching and extending similar sequencing fragment reads with the benchmark sequencing fragment reads to form a small contig.

[0027] Furthermore, step S7 specifically includes: there is thread overlap before and after in the generated thread identification table, which includes overlapping information of small contigs, constructing the before and after relationship of thread identifications through graph theory methods, constructing a multi-level and multi-node De Brujin Graph thread path, and performing sequence splicing according to the thread path to assemble the small contigs into contigs.

[0028] Furthermore, step S8 includes the following sub-steps:

[0029] S81: Extract the features of each contig base fragment;

[0030] S82: Compare the features of the contig base fragments, perform positional chimera on the corresponding contigs with the same features, and perform traversal error correction;

[0031] S83: Output the final assembly framework sequence scaffold.

[0032] A second-generation de novo assembly system based on gene numerical expression includes a memory management module, a sequencing fragment reads docking module, a small contig assembly module, a contig assembly module and a framework sequence scaffold assembly module, wherein the memory management module is used to realize base sequencing data management; the sequencing fragment reads docking module is used to parse the base sequence of the gene and customize the number, calculate the front and back centroids of two groups of specific lengths of each sequencing fragment read, build a B+ tree index through the base sequence and the number read ID of each sequencing fragment read, apply an artificial intelligence algorithm to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical fast arrangement algorithm to sort the centroid values ​​in sequence, and generate an ID docking table accordingly; the small contig assembly module is used to multi-threadedly assemble small contigs according to the ID docking table, and output a thread identification table, small contig fragments and their numbers and a single nucleotide variation (SNP) information table; the contig assembly module is used to construct a De Brujin Graph according to the thread identification and find a path to assemble small contigs The contigs are assembled into contigs; and the framework sequence scaffold assembly module is used to assemble framework sequence scaffolds based on the contigs.

[0033] The beneficial effects of the present invention are as follows: the present invention creatively proposes a second-generation de novo assembly method and system based on digital gene expression, which completes low-cost, high-performance assembly of second-generation sequencing. Compared with previous assembly algorithms, this method does not split reads into multiple kmer segments and construct a De Bruijn Graph through sequence overlap. Instead, it digitizes the sequencing bases through methods such as unsupervised machine learning, converting the order of overlap between reads into a numerical corresponding order. This not only improves the speed of read arrangement, but also avoids the generation of gap positions between kmer positions that may occur during the assembly process to a certain extent, reduces the number of times reads are repeatedly used for alignment, and greatly improves the calculation speed while ensuring accuracy. During the scaffold assembly process, instead of using the method of read re-alignment to increase computing power cost and IO overhead, a method of enzyme site identification-assisted assembly is adopted to extract enzyme site features of long fragment sequences, chimerize the gene files, and assemble them into longer fragments. During the assembly process, SNP variation information is also detected simultaneously. Using CUDA heterogeneous programming and adopting GPU acceleration strategies, CPU and GPU are coordinated to jointly calculate and improve parallel speed. In summary, this method is a second-generation assembly solution that combines high precision, low computing power requirements, and high efficiency. It significantly speeds up the bioinformatics analysis process for genome assembly analysis, helping bioinformatics analysts obtain accurate analysis results faster and in a shorter time.

[0034] Leveraging the rapid development of gene sequencing technology and bioinformatics, this method provides a powerful foundation for low-cost, rapid, and efficient genome assembly and analysis, truly realizing the storage, retrieval, and analysis of biological information using computer science and technology. Furthermore, bioinformatics is one of the core areas of natural science in the 21st century. As a key technology for DNA-seq analysis in bioinformatics, this method possesses significant technical advantages and market value in genomics, making it suitable for practical applications in genetic analysis and further technological advancement.

[0035] The key breakthroughs of this invention include: 1. It completely solves the memory-intensive problem of bioinformatics assembly. 2. It achieves complete traversal assembly, significantly reducing the risk of assembly errors. 3. It utilizes heterogeneous asynchronous computing to significantly improve computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying any creative work.

[0037] Figure 1 A flowchart for the invention;

[0038] Figure 2 It is a De Bruijn Graph path bubble structure graph;

[0039] Figure 3 It is a De Bruijn Graph path with a single in- and out-degree structure.

[0040] Figure 4 It is a De Bruijn Graph path ring structure graph;

[0041] Figure 5 This is a schematic diagram of head alignment;

[0042] Figure 6 This is a schematic diagram of tail alignment or staggered alignment;

[0043] Figure 7 Flowchart for small contig assembly;

[0044] Figure 8 Flowchart for contig assembly;

[0045] Figure 9 Assemble a flow chart for the framework sequence scaffold;

[0046] Figure 10 This is a system block diagram of the present invention. DETAILED DESCRIPTION

[0047] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions of the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations.

[0048] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings.

[0049] The following embodiments of the present invention are described in detail with reference to the accompanying drawings. In the absence of conflict, the following embodiments and features in the embodiments may be combined with each other.

[0050] Example 1:

[0051] See Figure 1 , a second-generation de novo assembly method based on gene digital expression, comprising the following steps:

[0052] S1: Base sequencing data management is achieved through soft disk array RAID;

[0053] S2: Gene analysis base sequence and custom numbering;

[0054] S3: Calculate the front and back centroid values ​​of two groups of specific lengths for each sequencing fragment read;

[0055] S4: Construct a B+ tree index based on the base sequence and read ID of each sequencing fragment read;

[0056] S5: Apply artificial intelligence algorithms to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical rapid permutation algorithm to sort the centroid values ​​in sequence, and generate an ID docking table based on this;

[0057] S6: Multi-threaded assembly of small contigs based on the ID docking table, outputting a thread identification table, small contig fragments and their numbers, and a single nucleotide variant (SNP) information table;

[0058] S7: Construct De Brujin Graph based on thread identifiers and find paths, assemble small contigs into contigs;

[0059] S8: Assemble scaffolds based on contigs.

[0060] In this embodiment, step S1 specifically includes: through a built-in RAID card or NFS protocol, the base sequencing data stored in the disk is transferred in a pseudo-virtual memory composed of a high IOPS hard disk in the form of a data stream, and is stored in a dispersed manner, wherein the base sequencing data includes centroid data and gene files. Step S2 specifically includes: reading the gene file from a high random read and write capability solid-state drive through the CPU to obtain multiple gene sequencing fragment reads, and numbering each gene sequencing fragment read to obtain a random number ID, and then transmitting the gene sequencing fragment reads and the random number ID to the GPU through the CUDA stream. Step S3 specifically includes: using the GPU to perform macro-layer identification on the base sequencing data, removing n bases before and after each sequencing fragment read (n value is large and small), obtaining two groups of four split front and back layer fragments, and using unsupervised machine learning to express the two groups of front and back layer sequences as centroid values ​​calculated by hard clustering, such as Kmeans clustering. The step S4 specifically includes: two sets of front and back centroid values ​​of specific lengths, constructing a B+ tree index with the base sequence of each sequencing fragment read and its corresponding numbered read ID to achieve random extraction of the base sequence. The step S5 specifically includes: after obtaining the front and back centroid values ​​of two sets of different fragments of all sequencing sequences, using an artificial intelligence algorithm, taking the two sets of front and back layer centroids and sequence lengths of the two sequences as input vectors, and whether the two sequences can match as the output vector, applying a deep learning neural network (DNN) to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, and generating a highly accurate ID docking table based on this, or using a numerical quick sort algorithm to sort the centroids in sequence, and generating an ID docking table based on this.

[0061] Furthermore, steps S1 to S5 can be regarded as the docking of sequencing fragment reads. Read docking can be roughly divided into three steps, the purpose of which is to obtain the read order that is preliminarily arranged according to the overlapping parts of the bases. The first step is to implement base sequencing data management, and transfer the data stored on the disk into the virtual memory through the built-in RAID card or NFS protocol, or use the virtual memory that comes with the Windows system. At the same time, the GPU is called to start the unsupervised machine learning kernel function, and each sequencing fragment read is marked with the front and back layer centroids respectively, or the front and back centroid values ​​of two groups of specific lengths of each sequencing fragment read are calculated. The second step uses the base sequence and numbered read ID of each sequencing fragment read to build a B+ tree index, and quickly sort the centroid values ​​in sequence. The third step involves pairing the randomly numbered IDs corresponding to the front- and back-end centroids with high similarity by setting a threshold or calculating data similarity. Alternatively, the two sets of front- and back-end centroids and sequence lengths of the two sequences are used as input vectors, and the output vector indicates whether the two sequences match. A deep learning neural network (DNN) is then applied to develop a matching model that aligns the front- and back-end centroids with the base sequence, resulting in a more accurate ID docking table. This provides the approximate order of all sequencing reads, i.e., the OLC (Overlap-Layout-Consensus) base sequence. For small genomic data, a multi-threaded approach can be used to first perform a full sequence I / O on each docking table, deleting unmatched read IDs and optimizing the ID docking table to obtain a completely correct clean ID docking table.

[0062] The docking of sequencing reads requires the use of a gene numerical expression method, which is mainly divided into two levels:

[0063] The first level is the macro-layer identification layer. The sequencing fastq files (gene files) after being downloaded from the machine are pre-processed by quality control to remove adapters and low-quality bases to obtain high-quality clean data. For the pre-processed fastq files (i.e., clean data), the method of breaking the sequencing fragment reads into fixed-length k-mers in the traditional DBG assembly algorithm is abandoned. The R1 file and R2 (gene file) files are respectively macro-layer identified, and m and n bases are removed from the front and back of a sequencing fragment read, respectively. Among them, the values ​​of m and n are different, and two groups of four split front and back layer fragments are obtained. Unsupervised machine learning is used to express the two groups of front and back layer sequences as the centroid values ​​calculated by Kmeans, and the front layer centroid and the back layer centroid are used to jointly identify a sequencing fragment read. For example, a preprocessed clean read is ATTTCGGGTCGAAA, 14 bp in length. Three bases are trimmed from the beginning and end, splitting it into ATTTCGGGTCG and TCGGGTCGAAA. The centroids of the front and back layers are calculated separately. The front and back layers can then collectively identify a single sequencing read. Specifically, in the macro-layer centroid calculation, A and T are designed as the positive and negative vertical axes, and G and C are designed as the positive and negative horizontal axes. Considering the position of the bases, each base in the split front and back sequences is represented as a point in a coordinate system, so that the base coordinates are scattered along the axes of the Cartesian coordinate system. For example, in the front layer sequence: TGCTT, T is the first base, with coordinates (0, -1), G is the second base, with coordinates (2, 0), and so on. The first front layer sequence can be converted into a coordinate array of {(0, -1), (2, 0), (-3, 0), (0, -4), (0, -5)}. Kmeans clustering was then used, using the mean squared error of the Euclidean distances between sample points as the criterion function. The centroid of this set of coordinate points was calculated, which served as the eigenvalue for the front-layer sequence. Similarly, the eigenvalue for the back-layer sequence was calculated. Weights were adjusted during the calculation process to ensure that the front-layer centroid was more sensitive to changes in the first 20 base pairs or so, while changes in a few bases in the remaining sequence did not affect the final front-layer centroid value. The back-layer centroid was more sensitive to changes in the last 20 base pairs or so, while changes in a few bases in the remaining sequence did not affect the final back-layer centroid value. The goal was to use the same centroid to represent as many similar sequencing reads as possible.

[0064] The method of extracting sequence feature identifiers of front and back layers of sequencing fragment reads of the present invention is not limited to the K-means clustering algorithm, but also includes other clustering algorithms, such as the maximum expectation clustering algorithm based on the Gaussian mixture model, the mean shift clustering algorithm, and the Fourier coefficient extraction feature, etc. The preferred embodiment of the present invention is to use the Kmeans clustering analysis algorithm to perform centroid calculation.

[0065] The second level is the centroid quick sort layer. After calculating the macro-level centroids, the resulting file includes the base sequences of the sequencing reads and the centroids of the preceding and following layers. These files are custom-numbered as read IDs. The centroids of the preceding and following layers corresponding to the read IDs are then numerically quick sorted by size to form a quick sort table. Using a built-in RAID card (disk array array card) or the NFS protocol (Network File System), the centroid data and gene files stored on disk are streamed into a pseudo-virtual memory consisting of high-IOPS (Input / Output Operating Per Second) hard drives. These data are then distributed and stored, forming a soft RAID. A B+ tree index is constructed between the read IDs of the fastq files and the base sequences corresponding to the randomly stored IDs in the soft RAID, facilitating the rapid retrieval of base sequences from the disk based on the read ID.

[0066] For large genomes, extraction is performed based on soft RAID and special hard drives. For small genome assembly, all intermediate files and sequencing data can be placed in the virtual memory of the Windows system for easy retrieval at any time.

[0067] See Figure 7 , the step S6 specifically includes: multi-threaded reading of the numbered read ID in the ID docking table, using the B+ tree index to extract the corresponding numbered read ID from the hard disk or virtual memory respectively, and matching them with the No. 1 sequencing fragment reads in sequence, matching and extending similar sequencing fragment reads with the No. 1 sequencing fragment reads to form a small contig. The purpose of assembling small contigs is to obtain highly accurate small contig fragments and SNP information tables. First, multi-threaded reading of the read ID in the ID docking table, using the B+ tree index to extract the corresponding base sequences from the hard disk or virtual memory respectively, and matching them in sequence, matching and extending extremely similar sequencing fragment reads to form a small contig and outputting the small contig's unique thread identifier, the assembled small contig, and the SNP information table.

[0068] See Figure 8, step S7 specifically includes: there is thread overlap before and after in the generated thread identification table, which includes overlapping information of small contigs, constructing the front and back relationship of the thread identification through graph theory methods, constructing a multi-level and multi-node De Brujn Graph thread path, and performing sequence splicing according to the thread path to assemble the small contigs into contigs. Further, the multi-threaded genomic contig assembly can be roughly divided into three steps, the purpose of which is to obtain high-accuracy small contig fragments and SNP information tables. First, the ID docking table with the complete thread identification and the scattered numbered reads ID in the thread identification are read and replaced with the thread identification in the docking table. Next, a De BruijnGraph directed graph is constructed according to the matching information of the thread identification, wherein each path contains the order of contig assembly. Finally, the small contigs corresponding to the thread identification are extracted according to the path in the De Bruijn Graph, assembled into contigs, and the information such as the variant base site and variant base type in the SNP docking table is modified.

[0069] Furthermore, steps S6 and S7 require the use of multi-threaded methods for constructing gene contigs and SNP variant detection. The multi-threaded method for constructing gene contigs can be generally divided into three levels: The first level is the ID docking level: After the centroid values ​​of the sequencing reads are initially aligned, multi-threaded reading is performed to read the front-layer centroids and back-layer centroids of reads 1-n (assuming read 1 is 0.000001, 0.000001; 2.359876, 5.164616). From the generated quick sort results, the front-layer centroids are found that are similar to the back-layer centroid of read 1 within a certain error range (for example, ±0.002). The IDs of the reads in this section are extracted to generate an ID docking table. These reads are 1-3, 9, 18, 66, 77, and 105.

[0070] In addition to directly setting a specified threshold range (i.e., error range) when performing ID matching, the present invention can also use a similarity measurement algorithm to calculate the similarity between centroid values ​​and match IDs corresponding to centroid values ​​with high similarity. Similarity calculation methods include calculating the Euclidean distance, cosine similarity, Manhattan distance, etc. between the centroid points of the front and back layers.

[0071] To optimize the ID docking table and ensure accurate and non-redundant matching between reads, a deep learning neural network was trained on the dataset, correlating the corresponding patterns of centroid values ​​with base sequence similarity. Specifically, the centroid values ​​were calculated for each read N bases before and after trimming. The first and last 60 bases of each read (reads less than 60 bases were calculated based on the entire base sequence) were then extracted and the centroid values ​​were calculated for each read, resulting in four sets of centroid values. KMP was used to perform a pairwise match of all base sequences. If the base sequences aligned (similar), a TRUE value was returned; if they did not, a FALSE value was returned. After obtaining the centroid coordinate values ​​of the paired matching sequences (4 groups of 8 centroid data in total) and the lengths of the two sequences, the ten components (the centroid coordinates of the first sequence with the first m bases removed (x1, y1), the centroid coordinates of the first sequence with the last m bases removed (x2, y2), the centroid coordinates of the first sequence with the first n bases removed (x3, y3), the centroid coordinates of the first sequence with the last n bases removed (x4, y4), the centroid coordinates of the second sequence with the first m bases removed (x5, y5), the centroid coordinates of the second sequence with the last m bases removed (x6, y6), the centroid coordinates of the second sequence with the first n bases removed (x7, y7), the centroid coordinates of the second sequence with the last n bases removed (x8, y48), the length of the first sequence Length1, the length of the second sequence Length2) are transformed to within [-1, 1] and used as the input vector x of the data set. The results of the matching calculation determine whether two sequences are similar or dissimilar. If similar, y = 1; if dissimilar, y = 0. y serves as the output vector of the dataset. Before DNN training, the network weights are initialized using a random distribution with a mean of 0. During training, the input vector x is fed into the input layer of the DNN network for forward propagation. After passing through several hidden layers, the forward output y' is obtained. The error between y' and the true value y is then calculated. Starting from the output layer, the error is backpropagated to the first hidden layer, adjusting the network's connection weights in a manner that minimizes the error. This process is repeated until the training error converges to a minimum, resulting in the desired DNN network weights.

[0072] After the DNN model training is completed, the target centroid file is identified to determine the one-to-many connection relationship between different sequences and output the ID docking table of the base sequence.

[0073] The second level is the sequence matching layer: multi-threaded reading of the numbered read IDs in the ID docking table, using the B+ tree index to extract the corresponding numbered read IDs (3, 9, 18, 66, 77, 105) from the hard disk or virtual memory, and matching them with the sequencing fragment read No. 1 in turn. The very similar sequencing fragment reads are matched and extended with No. 1 to form small contigs.

[0074] The specific multi-threaded implementation is as follows: When a thread begins calculation, a label is placed before the numbered read IDs in the ID docking table within the thread's scope, indicating a one-to-one match of base sequences 1-n. When the first thread performs calculations to match base sequence 1 with other docking base sequences, a label 1-1 is placed before the first row of the ID docking table (1-3, 9, 18, 66, 77, 105) as the thread identifier, indicating the first thread's second calculation. The corresponding sequences are then extracted based on the numbered read IDs using the B+ tree index and matched one-to-one. If the bases of sequence fragments 3 and 66 are largely identical to those of sequence 1, with only a few SNP variations, sequences 1, 3, and 66 are considered to overlap and can be merged into one sequence. If sequences 9, 18, and 105 have a significant number of base differences in their overlapping regions with sequence 1, these three sequences are discarded. When number 3 can match number 1, the ID docking table of sequence number 3 is continued to be searched. It is found that the third thread has started calculating the thread identification number 3-1 before the ID docking table. The calculation is stopped and the thread identification number {1-1: 3-1, 66, 77} of this calculation and the base sequence number 1 corresponding to the identification are output. The sequence of this small overlapping group samll contig is identified as contig1-1.

[0075] After completing the calculation for sequence read 1, the first thread starts a second calculation, searching for sequence reads not identified in the ID docking table, such as 9-101 and 102, and prefixing them with thread identifier 1-2. If reads 101 and 102 match read 9 and neither has been identified by another thread, the first thread extracts the corresponding base sequences, such as 101-1007 and 1008, and 102-1009 and 1010, and then splices them together. The calculation ends when a sequence with an identifier is found. If read 1010 has already been calculated by threads 8-15, the calculation ends and outputs the identifiers {1-2: 8-15, 1007, 1008, 1009} and the concatenated sequence of reads 9, 101, and 102. This small contig is identified as contig 1-2. The calculation continues in this order.

[0076] After all ID matching tables are calculated, all ID matching tables marked with thread identifiers are retrieved and the output thread identifier table is optimized. For example, if threads 4-7 and 8-9 have started the calculation of read 66 and read 77, respectively, the identifier table {1-1: 3-1, 66, 77} is corrected to {1-1: 3-1, 4-7, 8-9}. The thread identifier also represents the number of the output small contig.

[0077] The third layer is the graph-theoretic concatenation layer. The generated identifier table contains overlapping threads, including information about small contigs. Using graph-theoretic methods, the preceding and following relationships between thread identifiers are constructed, forming a multi-level, multi-node De Bruijn Graph thread path. Further sequence concatenation is performed along this path to form longer contig segments.

[0078] Specifically, after obtaining a series of thread identification tables, the De Bruijn Graph path diagram is constructed as follows: Figures 2 to 4 As shown. Among them, Figure 2 Bubble-like structure, Figure 3 It is a single-degree structure. Figure 4 For circular structures, for bubble-like graphs, they are assembled separately according to different paths. If the final assembly results of multiple paths are the same, they are merged and deduplicated and then output uniformly. If they are different, they are output separately. For graphs with a single degree of in-and-out, the assembly results are output separately according to different paths. For circular structures, this indicates that there may be repeated sequences. After assembly according to different paths, they are output separately once as repeated sequence files, which are focused on in the later assembly. When matching and assembling according to the ID docking table, the bases before and after the centroid value calculation are too long. There may be 10 to 20 bp before and after that can be matched, but the range represented by the centroid value is too large, so there is no docking problem in the ID docking table. For the assembly of this type of overlapping group contig fragments, it is necessary to generate centroid values ​​of different lengths. DNN (deep neural network) is used to calculate the matching relationship between two sets of centroid values ​​representing inconsistent base lengths, so as to optimize and generate a more accurate ID docking table.

[0079] For the multi-threaded assembly in the present invention, the above method is suitable for biological assembly contigs with large genomic data. If the genomic data is small, for the generated ID docking table, a multi-threaded solution can be first used to perform a full sequence IO on each docking table, delete the unmatched reads ID, optimize and delete the ID docking table to obtain a completely correct clean ID docking table, and then construct a DeBruijn Graph based on the entire ID docking table. Finally, the base sequence is extracted according to the path for splicing and assembly to generate contig fragments.

[0080] Furthermore, a method for detecting SNP variations: In the present invention, in addition to completing the de novo assembly of gene sequences, SNP (single nucleotide) variations of genes can also be detected. When constructing overlapping contigs in a multi-threaded manner, the SNP site information can be output by extracting sequencing fragment reads and matching them one by one: when assembling small overlapping contigs in a multi-threaded manner, the SNP site information is output for the first time, including the identifier, full length, SNP variation base site of the small overlapping contig, the base of the output overlapping contig at the variation site, the base replaced after the variation occurs, and other information. After the graph theory construction is completed, when extracting small overlapping contigs to assemble overlapping contigs, the SNP variation information must be adjusted accordingly.

[0081] Specifically, when extracting small contigs based on the path, the alignment is as follows: Figure 5 and Figure 6 Two cases are shown, where Figure 5 For head alignment, Figure 6 For tail alignment or staggered alignment, Figure 5 K1, K2, and K3 are the variant sites in the original SNP information table. If the alignment is head alignment, the original SNP variant information is directly output. If the variant sites of the reference sequence and the remaining sequences overlap, the replaced base information is merged and output. Figure 6 Where K1, K2, K3, K4, and K5 represent the variant sites in the original SNP information table, L represents the full length of the reference sequence, and O represents the overlap between the upper and lower bases after alignment. If the alignment is tail-aligned (staggered), the variant site in the reference sequence remains unchanged, while the variant sites in the remaining sequences are adjusted to L-O+K3, L-O+K4, and L-O+K5, respectively. If the variant sites in the reference sequence and the remaining sequences overlap, the replaced base information is merged and output.

[0082] The overall invention provides CPU and GPU computing environment detection functions for small gene assembly, and provides a CPU-adapted assembly method for personal computers without GPUs. The overall algorithm can use the CPU's AVX-512 deep learning instruction set to replace short fragments for more complex branch calculations.

[0083] See Figure 9 , the step S8 includes the following sub-steps:

[0084] S81: Extract the features of each contig base fragment;

[0085] S82: Compare the features of the base fragments of the contigs, perform positional chimera on the corresponding contigs, and perform traversal error correction;

[0086] S83: Output the final assembly framework sequence scaffold.

[0087] Furthermore, framework sequence scaffolds are assembled based on contig contigs. After obtaining contig sequences assembled from R1 and R2 (gene files), the bases of R2 are inverted to convert the sequence to a 5'-3' sequence. Specified enzyme digestion methods, such as 5'GGGCC^C 3' and 5'G^GATCC3', are used to identify sites. Single or double enzyme digestion can be used to tag all contigs, and a site tag interval table is output for each contig. Finally, contigs with the same site intervals are further assembled based on the number of intervals, forming longer framework sequence scaffolds. The SNP information table is then corrected in one step. In addition to using restriction digestion sites as contig features, deep learning methods such as bidirectional recurrent neural networks and transformer multi-head attention mechanisms can be used to generate a recurrent feature matrix for each contig. Contigs with similar features are then extracted for framework sequence scaffolding.

[0088] Step S8 of the present invention adopts a method of enzyme cutting site recognition-assisted scaffold assembly (framework sequence). This method discards all sequencing fragment reads and re-attaches them to the overlapping group contig after assembly, and then uses pair-end information to assemble the framework sequence scaffold. The method uses enzyme cutting site recognition-assisted assembly to identify different enzyme cutting sites and output the number of marker site intervals, thereby assembling the overlapping group contig into a framework sequence scaffold.

[0089] Specifically, in paired-end sequencing, the R1 file typically contains the 5'-3' end sequencing results, while the R2 file typically contains the 3'-5' end sequencing results. After obtaining the assembled contig sequences for R1 and R2, the bases in R2 are inverted to convert them to 5'-3' end sequences. All contigs are then site-tagged using designated enzyme recognition sites, and a site-tagged table for each contig is output. For example, if the base site cleaved by Apa I is 5'GGGCC^C3', the number of gaps between the 5'GGGCC^C 3' sequence in each contig is identified and output as a feature for each contig. During assembly, contigs with the same number of gaps are extracted and assembled into a scaffold. In addition to single enzyme digestion, two enzyme digestion sites can also be identified simultaneously, and the number of gaps between the two enzyme digestion sites can be expressed as positive or negative numbers for comparison.

[0090] The method for assembling scaffold sequences in the present invention mainly relies on the characteristics of contig sequences. In addition to using enzyme cleavage sites as features for assembly, a bidirectional recurrent neural network (Bi-RNN) or a transformer multi-head attention mechanism can also be used to extract features from contig sequences forward and backward to form a recurrent feature network, and contig assembly is performed based on similar features.

[0091] See Figure 10Based on the same inventive concept, the present invention also proposes a second-generation de novo assembly system based on gene numerical expression, which is used to implement the above-mentioned second-generation de novo assembly method based on gene numerical expression. The system includes a memory management module, a sequencing fragment reads docking module, a small contig assembly module, an overlapping contig assembly module and a framework sequence scaffold assembly module, wherein the memory management module is used to realize base sequencing data management; the sequencing fragment reads docking module is used to parse the base sequence of the gene and customize the numbering, calculate the front and back centroids of two groups of specific lengths of each sequencing fragment read, and construct a B+ tree index through the base sequence and number readID of each sequencing fragment read, apply an artificial intelligence algorithm to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical fast arrangement algorithm to sort the centroid values ​​in sequence, and generate an ID docking table accordingly; the small contig assembly module is used to assemble small contigs in multiple threads according to the ID docking table, and output a thread identification table, a small contig assembly module, and a small contig assembly module. contig fragments and their numbers and single nucleotide variation SNP information table; the contig assembly module is used to construct the De Brujin Graph and find the path according to the thread identifier, and assemble the small contigs into contigs; the framework sequence scaffold assembly module is used to assemble the framework sequence scaffold based on the contig contig.

[0092] As can be imagined, the purpose of the present invention is that genome assembly is a core issue in the field of bioinformatics and the starting point for a series of studies such as gene function analysis and downstream population genetic and structural differences. Based on a coordinate offset algorithm and an unsupervised machine learning hard clustering algorithm, the present invention uses a double-layer centroid to identify read sequences. Thread labels are used instead of kmer sequences split from sequencing files, and graph theory methods are used to construct Hamiltonian or Euler paths to find the order between reads. CUDA data streams are used in conjunction with asynchronous operations to avoid GPU (Graphics Processing Unit) waste and memory overhead, thereby improving computing speed. CUDA heterogeneous programming ensures high-speed communication between the CPU (Central Processing Unit / Processor) and GPU, enabling collaborative computing. This enables highly accurate and efficient de novo assembly of gene sequencing data, reducing the time cost of genome analysis. The CPU's AVX-512 instruction set is used to improve the efficiency of complex multi-branch computations.

[0093] The present invention has at least the following technical effects:

[0094] The present invention creatively proposes a second-generation de novo assembly method and system based on digital gene expression, achieving low-cost, high-performance assembly of second-generation sequencing. Compared with previous assembly algorithms, this method does not split reads into multiple kmer segments and construct a De Bruijn Graph based on sequence overlap. Instead, it uses unsupervised machine learning and other methods to digitize the sequencing bases and convert the order of overlap between reads into a numerical quick sort order. This not only improves the speed of read arrangement but also, to a certain extent, avoids the generation of gap positions between kmer positions during the assembly process, reduces the number of times reads are repeatedly used for alignment, and greatly improves the computing speed while ensuring accuracy. During the scaffold assembly process, instead of using the method of read re-alignment to increase computing power cost and I / O overhead, a restriction enzyme site-assisted assembly method is adopted to extract restriction enzyme site features of long fragment sequences, chimerize the gene files, and assemble them into longer fragments. During the assembly process, SNP variation information is also detected. Using CUDA heterogeneous programming and GPU acceleration strategies, CPU and GPU are coordinated to achieve parallel operation speed. In summary, this method is a second-generation assembly solution that combines high precision, low computing power requirements, and high efficiency. It significantly speeds up the bioinformatics analysis process for genome assembly analysis, helping bioinformatics analysts obtain accurate analysis results faster and in a shorter time.

[0095] Leveraging the rapid development of gene sequencing technology and bioinformatics, this method provides a powerful foundation for low-cost, rapid, and efficient genome assembly and analysis, truly realizing the storage, retrieval, and analysis of biological information using computer science and technology. Furthermore, bioinformatics is one of the core areas of natural science in the 21st century. As a key technology for DNA-seq analysis in bioinformatics, this method possesses significant technical advantages and market value in genomics, making it suitable for practical applications in genetic analysis and further technological advancement.

[0096] The key breakthroughs of this invention include: 1. It completely solves the memory-intensive problem of bioinformatics assembly. 2. It achieves complete traversal assembly, significantly reducing the risk of assembly errors. 3. It utilizes heterogeneous computing to significantly improve computational efficiency.

[0097] It should be noted that for the aforementioned embodiments, for simplicity of description, they are all expressed as a series of action combinations, but those skilled in the art should be aware that this application is not limited by the order of the actions described, because according to this application, some steps can be performed in other orders or simultaneously. Secondly, those skilled in the art should also be aware that the embodiments described in this specification are preferred embodiments, and the actions involved are not necessarily required by this application.

[0098] The above embodiments describe the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely illustrative of the principles of the present invention. Without departing from the spirit and scope of the present invention, modifications and variations made by those skilled in the art without departing from the spirit and scope of the present invention should be within the scope of protection of the appended claims.

Claims

1. A second-generation de novo assembly method based on digital gene expression, characterized in that: The steps include: S1: Base sequencing data management is achieved through soft disk array RAID; S2: Gene analysis base sequence and custom numbering; S3: Calculate the front and back centroid values ​​of two groups of specific lengths for each sequencing fragment read; S4: Construct a B+ tree index based on the base sequence and read ID of each sequencing fragment read; S5: Apply artificial intelligence algorithms to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical rapid permutation algorithm to sort the centroid values ​​in sequence, and generate and optimize the ID docking table based on this; S6: Multi-threaded assembly of small contigs based on the ID docking table, outputting a thread identification table, small contig fragments and their numbers, and a single nucleotide variant (SNP) information table; S7: Construct De Brujin Graph based on thread identifiers and find paths, assemble small contigs into contigs; S8: Assemble scaffolds based on contigs.

2. A second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S1 specifically includes: using a built-in RAID card or NFS protocol, the base sequencing data stored in the disk is transferred in a pseudo-virtual memory composed of a high-IOPS hard disk in the form of a data stream, and stored in a dispersed manner, wherein the base sequencing data includes centroid data and gene files.

3. A second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S2 specifically includes: reading the gene file from the solid-state drive with high random read and write capabilities through the CPU to obtain multiple gene sequencing fragment reads, and numbering each gene sequencing fragment read to obtain a random number ID, and then transmitting the gene sequencing fragment reads and random number ID to the GPU through the CUDA stream.

4. A second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S3 specifically includes: using the GPU to perform macro-layer identification on the base sequencing data, removing m and n bases before and after each sequencing fragment read, respectively, where the values ​​of m and n are different, to obtain two groups of four split front and back layer fragments, and using unsupervised machine learning to express the two groups of front and back layer sequences as centroid values ​​calculated by hard clustering.

5. The second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S4 specifically includes: after calculating the front and back centroid values ​​of two groups of specific lengths, constructing a B+ tree index for the base sequence of each sequencing fragment read and its corresponding number read ID to achieve random extraction of the base sequence.

6. A second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S5 specifically includes: after obtaining the front and back centroid values ​​of two groups of different fragments of all sequencing sequences, using an artificial intelligence algorithm, taking the two groups of front and back layer centroids and sequence lengths of the two sequences as input vectors, and whether the two sequences can match as the output vector, applying a deep learning neural network DNN to obtain a matching model of the relationship between the front and back centroid values ​​and the base sequence alignment, and generating a highly accurate ID docking table based on this, or using a numerical quick sort algorithm to sort the centroids in sequence, and generating an ID docking table based on this.

7. A second-generation de novo assembly method based on gene digital expression according to claim 6, characterized in that: The artificial intelligence algorithm of step S5 includes the following sub-steps: S51: transform the two sets of front and back centroid values ​​of the two sequences and their sequence lengths to within [-1, 1] as the input vector x of the data set; S52: Use the KMP algorithm to align the two sequences. If they can be aligned, it returns TRUE, y=1; if they cannot be aligned, it returns FALSE, y=0, and y is used as the output vector of the data set. S53: Before training, first set a random distribution with a mean of 0 to initialize the network weights. During training, the input vector x is first sent to the input layer of the DNN network for forward propagation. After passing through several hidden layers, the forward output y' is obtained; S54: Calculate the error between y' and the true value y, and then backpropagate the error from the output layer to the first hidden layer, adjusting the connection weights of the network in the direction of reducing the error; S55: The above weight adjustment process is continuously repeated until the training error converges to a minimum value, that is, the DNN network weights that meet the requirements are obtained.

8. The second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S6 specifically includes: multi-threaded reading of the numbered read IDs in the ID docking table, using a B+ tree index to extract the corresponding numbered read IDs from the hard disk or virtual memory, and matching them with the benchmark sequencing fragment reads in sequence, matching and extending similar sequencing fragment reads with the benchmark sequencing fragment reads to form small contigs.

9. The second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S7 specifically includes: there is thread overlap before and after in the generated thread identification table, which includes overlapping information of small contigs, constructing the before and after relationship of thread identifications through graph theory methods, constructing a multi-level and multi-node De Brujin Graph thread path, and performing sequence splicing according to the thread path to assemble the small contigs into contigs.

10. A second-generation de novo assembly method based on gene digital expression according to claim 1, characterized in that: The step S8 includes the following sub-steps: S81: Extract the features of each contig base fragment; S82: Compare the features of the contig base fragments, perform positional chimera on the corresponding contigs with the same features, and perform traversal error correction; S83: Output the final assembly framework sequence scaffold.

11. A second-generation de novo assembly system based on digital gene expression, for implementing the second-generation de novo assembly method based on digital gene expression according to any one of claims 1 to 10, characterized in that: The invention comprises a memory management module, a sequencing fragment reads docking module, a small contig assembly module, a contig assembly module and a scaffold assembly module, wherein the memory management module is used to realize base sequencing data management; the sequencing fragment reads docking module is used to parse the base sequence of the gene and customize the numbering, calculate the front and back centroids of two groups of reads of specific lengths for each sequencing fragment, build a B+ tree index through the base sequence and number read ID of each sequencing fragment read, apply an artificial intelligence algorithm to obtain a matching model of the alignment relationship between the front and back centroid values ​​and the base sequence, or use a numerical fast permutation algorithm to sort the centroid values ​​in sequence, and generate an ID docking table accordingly; the small contig assembly module is used to multi-thread assemble small contigs according to the ID docking table, and output a thread identification table, small contig fragments and their numbers and a single nucleotide variation SNP information table; the contig assembly module is used to construct a De Brujin Graph according to the thread identification and find a path, and assemble small contigs into a scaffold. The contigs are assembled into contigs; and the framework sequence scaffold assembly module is used to assemble framework sequence scaffolds based on the contigs.

Images