Methods, systems, devices, and media for assembling third-generation sequencing data based on clustering and graph construction.

Abstract

This application discloses a method, system, device, and medium for assembling third-generation sequencing data based on clustering and graph construction, belonging to the field of biological sequence processing technology. The method includes the following steps: obtaining sequence similarity based on the third-generation sequencing data; clustering the sequences using a clustering algorithm to obtain different clusters; assembling sequences within the same cluster using a graph construction method to obtain intra-group consensus sequences; merging all intra-group consensus sequences from different clusters and assembling them again using a graph construction method to obtain inter-group consensus sequences. The technical solution of this application effectively reduces the complexity of assembly through a "sequence clustering-intra-group assembly-inter-group assembly" strategy, and ensures the accuracy and integrity of the consensus sequences through the optimal overlap graph algorithm and deep pruning, exhibiting significant technical advantages.

Description

Technical Field

[0001] This application relates to the field of biological sequence processing technology, and in particular to a method, system, device and medium for assembling third-generation sequencing data based on clustering and graph construction. Background Technology

[0002] Amplicon sequencing is a key technique in molecular biology research. It specifically amplifies target gene regions using PCR (polymerase chain reaction), allowing for sequencing analysis of those regions. It is commonly used to reveal the diversity of target genes in a sample; for example, amplification and detection of 16S / 18S / ITS rRNA genes can reveal the structure of environmental microbial communities. Similarly, in clinical research, it can sequence specific genes in human samples (such as cancer tissue) to reveal disease-related gene variations.

[0003] With the maturation of third-generation sequencing technologies such as PacBio HiFi and Oxford Nanopore, those skilled in the art can obtain thousands to tens of thousands of complete amplicon reads, aiming to reconstruct high-fidelity consensus sequences for every gene variant present in a sample from these raw sequencing reads. However, the bioinformatics field currently lacks high-fidelity assembly tools specifically designed for this type of sequencing data (i.e., long reads, high depth, high similarity, and extremely high length variability). Those skilled in the art typically only have access to two existing technical approaches not designed for this purpose: The first technical approach employs tools designed for de novo genome assembly, such as Canu or Flye. The core algorithms of these tools (e.g., overlap-layout-consensus or A-Bruijn graphs) are optimized for processing random, low-depth genomic data, designed to address repetitive sequences and extend contigs. However, when faced with amplicon (high-depth, highly similar) data, these algorithms tend to incorrectly collapse similar variants (e.g., different alleles) into a single or chimeric consensus sequence. Furthermore, the filtering mechanisms of these tools are highly likely to incorrectly exclude overly long (carrying large insertions) or overly short (carrying large deletions) valid variants in cancer, leading to the loss of crucial variant information.

[0004] The second technical approach is to use traditional multiple sequence alignment tools (such as MUSCLE or MAFFT). This method performs a global alignment of all reads and then extracts a consensus sequence. This approach can achieve good results when the sample is clean and the total number of reads is small, but the computational complexity of its algorithm increases exponentially or exponentially with the number of sequences, making it extremely time-consuming to process the high-throughput reads generated by third-generation sequencing. Furthermore, when dealing with sequences of extremely varying lengths, multiple sequence alignment tools can introduce a large number of artificial gaps, leading to alignment failures or the generation of meaningless consensus sequences.

[0005] The existing technical approaches mentioned above suffer from drawbacks such as low assembly efficiency, high chimerism rate, and inability to construct accurate consensus sequences when dealing with the assembly of large-scale third-generation amplicon data due to misaligned tool design and lack of scalability of algorithms. Therefore, there is an urgent need for an efficient assembly method specifically designed for this scenario. Summary of the Invention

[0006] To solve at least one of the above-mentioned technical problems, the technical solution adopted in this application is as follows.

[0007] The first aspect of this application provides a method for assembling third-generation sequencing data based on clustering and graph construction, comprising the following steps: Sequence clustering: Based on the third-generation sequencing data, sequence similarity is obtained, and clustering algorithms are used to cluster the sequences to obtain different clusters; Intra-group assembly: Sequences within the same cluster are assembled using a graph-based construction method to obtain intra-group consensus sequences; Inter-group assembly: All intra-group consensus sequences from different clusters are merged and assembled using the same graph construction method to obtain inter-group consensus sequences.

[0008] In some embodiments of this application, the following steps are further included before the sequence clustering step: Sequence filtering: Filter out sequences whose length is greater than a first preset threshold and whose length is less than a second preset threshold. The first preset threshold and the second preset threshold are determined based on the range of all sequence lengths in the third-generation sequencing data.

[0009] In some embodiments of this application, the mean and standard deviation of all sequence lengths in the third-generation sequencing data are first calculated, and the first preset threshold is set with reference to the following formula: Mean + (Standard deviation × Factor); The second preset threshold is set with reference to the following formula result: Mean - (Standard deviation × Multiple factor) The multiplier factor is 1 to 2.

[0010] In this application, the reference refers to a value directly obtained from the calculation result of the formula, or in the form of rounding, taking a similar value, etc.

[0011] In some embodiments of this application, the third-generation sequencing data is obtained by sequencing amplicon sequences using a third-generation sequencer, and the sequence length in the third-generation sequencing data has a expected length. The first preset threshold is set with reference to the following formula result: Expected length × lower limit of length range ratio; The second preset threshold is set with reference to the following formula result: Expected length × lower limit of length range ratio The lower limit of the length range ratio is 0.4 to 0.8, and the upper limit of the length range ratio is 1.5 to 2.5.

[0012] Amplicon sequencing is a highly targeted method used to analyze gene variations in specific genomic regions. Amplicon sequencing mainly includes 16S rRNA gene sequencing, 18S rRNA gene sequencing, ITS gene sequencing, and target region amplicon sequencing. Since amplicon sequencing sequences PCR amplification products, the expected sequencing sequence length can be calculated based on primer and target / target region information.

[0013] In some embodiments of this application, the sequence filtering step further includes filtering out sequences with a length less than a third preset threshold. In the previous step, the first preset threshold may be set too low, resulting in many short sequences remaining unfiltered. Therefore, a global minimum filtering length is further set. In some embodiments of this application, the third preset threshold is 100 bp, which also means that the assembly method of this application can utilize sequences longer than 100 bp.

[0014] In this application, those skilled in the art can employ any reliable sequence similarity calculation method, such as sequence sketching methods or k-mer methods. Sequence sketching methods include, for example, Strobemer (Randstrobes); k-mer methods include, for example, MinHash or FracMinHash, used to quickly calculate Jaccard similarity between sequences; another example is full k-mer profiles, which use complete k-mer frequency vectors to compare sequences by calculating cosine similarity or Jaccard similarity.

[0015] In some embodiments of this application, during the sequence clustering step, The sequence similarity is calculated based on Strobemer: a Randstrobes feature vector is generated for each sequence; an index of Strobemer features is constructed, and candidate sequence pairs sharing Strobemer are found; bidirectional similarity is calculated only for the selected candidate pairs.

[0016] In some specific embodiments of this application, in the step of generating Randstrobes feature vectors for each sequence, the strobe length k takes a value of 18 to 21, the minimum window offset w_min takes a value of 10 to 15, and the maximum window offset w_max takes a value of 30 to 40. In some preferred embodiments of this application, k=20, w_min=12, and w_max=35.

[0017] In some specific embodiments of this application, strobemers that occur too frequently are filtered out, and the filtering threshold can be based on a frequency ratio (e.g., 0.03 to 0.1) or an absolute number (e.g., 150). Only candidate pairs with a shared strobemer count greater than or equal to a minimum threshold (e.g., 3) are retained. In some preferred embodiments of this application, the frequency ratio threshold for high-frequency hashes is set to 0.05, the absolute threshold is set to 150, and the minimum threshold for the number of shared strobemers is set to 3.

[0018] In some embodiments of this application, the similarity threshold is calculated to be 0.05 to 0.1, for example, 0.08.

[0019] Furthermore, those skilled in the art can employ various clustering methods. Progressive density clustering is an effective method. Those skilled in the art can also use other clustering algorithms, such as DBSCAN, a classic density clustering algorithm; K-Means clustering, which can cluster the feature vectors; and hierarchical clustering, used to construct a hierarchical structure of clusters.

[0020] In this application, progressive density clustering is adopted, and the similarity threshold is gradually lowered at different stages. At each stage, core points and boundary points are identified for each cluster.

[0021] In some embodiments of this application, the minimum cluster size and core cluster size are dynamically calculated based on the total number of sequences to be clustered.

[0022] In some embodiments of this application, when performing multi-stage clustering, the initial threshold is relatively high (e.g., the 80th to 85th percentile of the similarity score), and the threshold is gradually reduced (e.g., step number = 3) until the final threshold is reached (e.g., the 60th to 70th percentile of the similarity score).

[0023] Furthermore, after staged clustering, unassigned sequences are attempted to be assigned to existing clusters or to form new clusters based on a lower similarity threshold.

[0024] In some embodiments of this application, the maximum number of sequences per cluster is limited (e.g., 200 to 600); while clusters that are too small (e.g., minimum size = 3) are filtered out.

[0025] In some embodiments of this application, intra-group assembly and / or inter-group assembly are performed based on the Best Overlap Graph (BOG). The BOG is an efficient graph simplification algorithm. Those skilled in the art can also use other graph construction algorithms within each "purified" cluster, such as: the standard Overlap-Layout-Consensus (OLC) graph: constructing a complete overlap graph without BOG simplification; or: the Debugging Graph (DBG): since the number of sequences within each cluster is greatly reduced, it becomes possible to construct small Debugging graphs within the cluster, which can be used to resolve more fine-grained mutations.

[0026] In some specific embodiments of this application, the construction of the overlap graph includes the following three stages: Overlapping phase: In this stage, the algorithm searches for overlapping regions between all sequences. Overlap is defined as the length of the sequence ends matching the beginning of another sequence. By comparing sequences, pairwise overlaps are identified. All found overlaps are recorded, generating an overlap graph where each node represents a sequence and edges represent the overlap between two sequences.

[0027] In some embodiments of this application, minimap2 is used to find overlapping regions. Minimap2 is currently an efficient overlap detection tool, but any tool that can generate similar PAF format overlap information can be used as an alternative.

[0028] Composition stage: After the overlap graph is constructed, the algorithm processes it to find an optimal path to arrange the sequences, reconstructing the original genome sequence as closely as possible. Typically, this step involves simplifying the overlap graph, removing unnecessary edges (such as redundant overlaps in repetitive regions), and identifying the optimal path. The simplified graph structure forms a "layout graph," which represents the sequence order on the genome.

[0029] In some embodiments of this application, the overlap graph construction steps are as follows: Preprocessing PAF: Identifying and filtering suspicious overlaps arising from multiple matches. Classifying overlap relationships and calculating the error rate; Create a node: Create a graph node for each sequence; Create edges: Create directed edges based on overlapping records, recording the source node, target node, overlap length, consistency, overlap type, etc. Assign edge endpoints: Determine the positions of the 5' and 3' ends of each edge on the source and target nodes; Handling inclusion relationships: Identify inclusion relationships between sequences, cluster the included sequences, and select the longest sequence as the representative sequence to simplify the graph structure; Assigning optimal edges: For each node, assign the optimal incoming and outgoing edges to its 5' and 3' ends respectively based on overlap consistency.

[0030] Identify mutually beneficial edges: Identify edge pairs that are the best of each other.

[0031] Identify the starting point of the path: Prioritize strict starting points (i.e., nodes without 5' entry edges but with 3' exit edges). If no strict starting point exists, use lenient starting points (i.e., active nodes with 3' edges that are neither contained nor chimeras). Path expansion: Starting from the origin, extend the path forward along the 3' direction edge of the reciprocity. During the expansion, skip nodes that have been marked as containing or nested, and avoid reusing nodes. Unitig Construction: Construct unitigs from the solved paths, merge the sequences represented by all nodes on the path according to their overlap, and for the remaining nodes, identify the top-level container nodes, contained nodes, and independent sequences, and construct container-type unitigs and independent unitigs respectively, filtering out unitigs with shorter lengths.

[0032] Consensus Phase: After sequencing, the final step is to generate a consensus sequence. The consensus sequence is determined based on the sequence information of overlapping regions and is used to correct errors in sequencing. By analyzing the overlapping regions of multiple reads, the consensus sequence for each base is determined, resulting in a more accurate genome sequence. Common consensus algorithms count the frequency of bases at each position and select the most frequent base as the consensus base for that position, thus generating the final genome sequence.

[0033] In some embodiments of this application, the algorithms or software used to achieve optimal overlap graph construction include, but are not limited to: Overlap Layout Consensus (OLC), Canu, WTDBG2, and Falcon.

[0034] In some embodiments of this application, after the intra-group assembly and / or inter-group assembly steps, a step of trimming and polishing the sequence is also included: Sequence alignment and depth calculation: The sequence before clustering is used as a reference sequence and compared with the consensus sequence within or between groups, and the depth of each base position is calculated. Sequence pruning: Prune regions with a depth lower than the fourth preset threshold; Sequence polishing: Correcting the bases at each position based on alignment quality fraction, depth, and / or number of bases.

[0035] In some embodiments of this application, the sequence pruning includes two pruning methods: end-to-end pruning and middle pruning, wherein: Trimming at both ends: Starting from both ends of the sequence, examine inwards and remove regions with a depth lower than "average depth × first depth ratio threshold". Preferably, the first depth ratio threshold is 0.1~0.2.

[0036] Intermediate pruning: Processing regions in the middle of the sequence whose depth is lower than "average depth × second depth ratio threshold". Preferably, the first depth ratio threshold is 0.3~0.5. The processing is selected from one of the following three modes: Marking mode: Marks only low-depth regions in the output file without changing the sequence; Segmentation pattern: The sequence is split into two independent sequences at the midpoint of the low-depth region; Deletion mode: Delete the low-depth region from the sequence.

[0037] In some embodiments of this application, the trimming and polishing steps include multiple rounds, with the sequence modified and polished in the previous round serving as the reference sequence for the next round of comparison.

[0038] It is worth noting that, as mentioned above, the technical solution of this application, especially the trimming and polishing steps, aims to improve the base fidelity of the consensus sequence. However, third-generation sequencing platforms (such as Oxford Nanopore) exhibit a systematic shortening bias when processing homopolymers with long continuous lengths (e.g., exceeding 7 bp). (For example, a true 8 bp homopolymer may be systematically reported as 7 bp.) This error is difficult to completely correct using majority voting or weighted voting methods.

[0039] To address this issue, this application may also include an optional, additional "homogeneous correction step," which, following the trimming and polishing steps, is used to re-correct the homopolymer regions in the polished consensus sequence. Possible implementations include, but are not limited to: Signal-based correction: For data that supports the original signal (such as ONT's FAST5 / POD5 format), the original electrical or current signal data corresponding to the homopolymer regions in the consensus sequence output from the trimming and polishing steps can be reread and analyzed. By analyzing the duration or morphology of the signal, the exact length of the homopolymer can be re-estimated at higher resolution, thereby correcting the sequence.

[0040] Statistical model-based correction: A homopolymer error model for a specific sequencing platform and base type (e.g., A / T / C / G) is pre-built using a training set (e.g., data from both third-generation sequencing and the gold standard Sanger sequencing). (For example, an 8bp 'A' homopolymer has a 90% probability of being measured as 7bp and a 10% probability of being measured as 8bp). During correction, this model is used to probabilistically correct the homopolymer length in the consensus sequence.

[0041] Correction based on hybrid sequencing: In applications that support hybrid sequencing data, short-read sequencing data (e.g., Illumina data) with high precision and low homopolymer errors are accepted. Final error correction of the long-read consensus sequence is performed by aligning the short-read data back to the long-read consensus sequence output from the trimming and polishing steps, utilizing the explicit readings of the short-read data in the homopolymer regions.

[0042] A second aspect of this application provides a system for assembling third-generation sequencing data based on clustering and graph construction, comprising the following modules: The data input module is used to: obtain the third-generation sequencing data; The sequence clustering module, connected to the data input module, is used to: obtain sequence similarity based on the third-generation sequencing data, and cluster the sequences using a clustering algorithm to obtain different clusters; The intra-group assembly module, connected to the sequence clustering module, is used to: assemble sequences within the same cluster using a graph-based construction method to obtain intra-group consensus sequences; The inter-group assembly module, connected to the intra-group assembly module, is used to: merge all intra-group consensus sequences from different clusters and assemble them using the same graph construction method to obtain inter-group consensus sequences.

[0043] In some embodiments of this application, the system further includes a sequence trimming and polishing module connected to the data input module, and also connected to the intra-group assembly module and / or inter-group assembly module, for implementing the following steps: Sequence alignment and depth calculation: The sequence before clustering is used as a reference sequence and compared with the consensus sequence within or between groups, and the depth of each base position is calculated. Sequence trimming: Trim regions with a depth lower than the fourth preset threshold; Sequence polishing: Correcting the bases at each position based on alignment quality fraction, depth, and / or number of bases.

[0044] A third aspect of this application provides an electronic device including at least one processor and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform any of the methods described in the first aspect of this application.

[0045] A fourth aspect of this application provides a non-transitory computer-readable storage medium storing computer instructions for causing the computer to perform any of the methods described in the first aspect of this application.

[0046] Compared with the prior art, this application has the following advantages: The method in this application pre-separates large-scale raw reads into multiple small, independent clusters by employing efficient Strobemer similarity calculation and progressive density clustering. Then, the optimal overlap graph algorithm is used to perform assembly independently within each "purified" cluster, and finally, all raw reads are used for alignment depth-based pruning and polishing.

[0047] Compared with existing technologies (such as Canu / Flye or MSA tools), the technical solution provided in this application not only significantly improves the computing speed and solves the scalability problem of high-throughput data (such as tens of thousands of reads) through this "divide and conquer strategy" (solving the infeasibility of MSA computation); but also, through the "separate first, then assemble" mode, it fundamentally avoids the erroneous "folding" of different variants (such as alleles), significantly reducing the chimerism rate (solving the assembly error defect of whole genome tools); finally, combined with flexible filtering and polishing, it ensures the preservation of key variant information (such as large fragment insertions / deletions in cancer) and extremely high base fidelity (such as the 99.99% consistency verified in Example 1), thereby significantly improving the overall accuracy of amplicon assembly.

[0048] The technical solution of this application can be used not only for detecting microbial communities and cancer gene mutations, but also for the following technical fields: 1. Quasispecies analysis: used to isolate and reconstruct multiple highly similar but key site-different viral strains (such as HIV, HCV, and COVID-19 variants) from a single patient sample with high fidelity.

[0049] 2. Immune Repertoire Analysis: Used to assemble highly variable antibodies or T-cell receptors (TCRs) from the V(D)J recombination region to accurately reconstruct the clonal type.

[0050] 3. Specific gene family analysis in metagenomics: used to specifically isolate and assemble a class of highly conserved but variable functional genes (such as antibiotic resistance genes, rpoB genes, etc.) from complex metagenomic samples.

[0051] It should be understood that the description in this section is not intended to identify key or essential features of the embodiments of this application, nor is it intended to limit the scope of this application. Other features of this application will become readily apparent from the following description. Attached Figure Description

[0052] The above and other objects, features, and advantages of exemplary embodiments of this application will become readily apparent from the following detailed description taken in conjunction with the accompanying drawings. Several embodiments of this application are illustrated in the drawings by way of example and not limitation, in which: Figure 1 This paper shows a length distribution diagram of the sequences used for assembly in Embodiment 1 of this application; Figure 2 This illustration shows a flowchart of the assembly of amplicon third-generation sequencing data based on clustering and graph construction in Embodiment 1 of this application; Figure 3 This illustrates the clusters formed by progressive density clustering in Embodiment 1 of this application; Figure 4 This illustration shows a schematic diagram of a system for assembling amplicon third-generation sequencing data based on clustering and graph construction in Embodiment 3 of this application; Figure 5 This illustration shows an improved system for assembling amplicon third-generation sequencing data based on clustering and graph construction, as described in Embodiment 4 of this application. Detailed Implementation

[0053] Unless otherwise stated, implied from the context, or as is customary in the art, all parts and percentages in this application are based on weight, and all testing and characterization methods used are concurrent with the filing date of this application. Where applicable, any patent, patent application, or disclosure relating to this application is incorporated herein by reference in its entirety, and its equivalent patent families are also incorporated herein by reference, particularly the definitions of relevant terms in the art disclosed in such documents. If any definition of a specific term disclosed in the prior art is inconsistent with any definition provided in this application, the definition provided in this application shall prevail.

[0054] To make the technical problems, technical solutions and beneficial effects solved by this application clearer, the following detailed description is provided in conjunction with embodiments.

[0055] The following examples are used to illustrate preferred embodiments of this application. Those skilled in the art will understand that the techniques disclosed in the examples represent technologies discovered by the inventors that can be used to implement this application, and therefore can be considered preferred embodiments of this application. However, those skilled in the art should understand from this specification that many modifications can be made to the specific embodiments disclosed herein, still yielding the same or similar results, without departing from the spirit or scope of this application.

[0056] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains, and all materials cited herein and referenced by them are incorporated herein by reference.

[0057] Those skilled in the art will recognize, or can learn through routine experimentation, many equivalents of the specific embodiments of the invention described herein. These equivalents will be included in the claims.

[0058] Unless otherwise specified, the experimental methods used in the following examples are conventional methods. Unless otherwise specified, the instruments and equipment used in the following examples are all conventional laboratory instruments and equipment; unless otherwise specified, the experimental materials used in the following examples were all purchased from conventional biochemical reagent stores.

[0059] Example 1: A method for assembling amplicon third-generation sequencing data based on clustering and graph construction. This embodiment is based on the assembly of a long-read amplicon sequencing dataset obtained using the Qitan Tech QNome-3841 nanopore gene sequencer. This dataset contains approximately 2209 ONT reads (sequencing sequences). The amplicon was obtained by amplifying the cDNA of the mouse lactoferrin gene incorporating iCre, IRES, and EGFP, with an expected length of approximately 2442 bp. The actual length distribution is shown in the figure below. Figure 1 As shown. (Refer to...) Figure 2 The flowchart shown mainly includes steps S1 to S6, which are explained in detail below: S1, length filtering of amplicon sequencing sequences This step filters amplicon sequencing sequences to remove sequences with abnormal lengths, ensuring data quality for subsequent analyses.

[0060] First, the mean length (941 bp in this example) and standard deviation (590 bp in this example) of the 2209 read segments were calculated.

[0061] Next, select the filter mode. There are two filter modes to choose from: (1) Specify length mode For amplicon sequencing, an expected length can be obtained based on the amplification primers and target. The minimum length threshold for the filtering range is set as: expected length × lower limit of the length range ratio (e.g., 0.6); the maximum length threshold for the filtering range is set as: expected length × upper limit of the length range ratio (e.g., 2.0).

[0062] (2) Adaptive mode The minimum length threshold for the filtering range is set to: mean - (standard deviation × factor) (e.g., 1.5 times); the maximum length threshold for the filtering range is set to: mean + (standard deviation × factor) (e.g., 1.5 times).

[0063] In this embodiment, the "specified length mode" is selected for filtering. The expected length is set to 2442bp, with a range of 0.6 to 2.0 times. The calculated filtering range is: [expected length (2442bp) × 0.6, expected length (2442bp) × 2.0], that is: [1465bp, 4884bp].

[0064] Meanwhile, the full length of the third-generation sequencing adapter is 53bp~75bp. To prevent the threshold from being too low and causing dimers to affect the assembly results, the threshold needs to be set to approximately twice the adapter length. To ensure successful assembly, this embodiment adopts an aggressive strategy, setting the minimum threshold to approximately twice the lower limit of the adapter length, 253×2=106 (bp), and finally setting the global minimum length lower limit to 100bp (even so, since 1465>100, this lower limit has no effect in this embodiment).

[0065] All sequences whose lengths are outside the filtering range ([1465bp, 4884bp]) are filtered out. In this embodiment, 340 sequences remain after filtering, with the length distribution as follows: Figure 1 The area defined by the two red vertical dashed lines.

[0066] S2, Sequence clustering based on Strobemer similarity and clustering. In this step, based on Strobemer (specifically implemented as Randstrobes) similarity calculation and progressive density clustering algorithm, highly similar sequences are clustered into the same cluster, preparing for subsequent intra-group assembly.

[0067] I. Three-stage Strobemer similarity calculation 1. Feature Generation Randstrobes feature vectors are generated from the 340 sequences obtained in step S1. In this embodiment, the strobe length k is set to 20, the minimum window offset w_min is set to 12, and the maximum window offset w_max is set to 35.

[0068] 2. Candidate Screening Build an index of strobemer features and find candidate sequence pairs that share strobemers. To improve efficiency, filter out high-frequency hashes (i.e., strobemers that appear too frequently), with the filtering threshold based on frequency ratio (e.g., 0.03 to 0.1) or absolute number (e.g., 150). Only retain candidate pairs that share a minimum threshold (e.g., 3) of strobemers.

[0069] In this embodiment, the frequency ratio threshold for high-frequency hashing is set to 0.05, the absolute threshold is set to 150, and the minimum threshold for the number of shared strobemers is set to 3.

[0070] 3. Similarity Calculation For selected candidate pairs only, calculate bidirectional containment similarity (i.e., the proportion of features contained in sequence B by sequence A, and the proportion of features contained in sequence A by sequence B). A threshold for similarity calculation is used (e.g., 0.05 to 0.1). The calculated similarity matrix is ​​stored in sparse matrix format to save memory.

[0071] In this embodiment, the similarity threshold is set to 0.08. At this threshold, the vast majority of potential homologous sequences can be recalled, while random similar sequences caused by sequencing errors are filtered out.

[0072] II. Progressive density clustering 1. Dynamic parameter calculation Based on the total number of filtered sequences, the minimum cluster size and core cluster size are dynamically calculated. Specifically, the process is as follows: (1) Select a step-decreasing ratio coefficient based on the total number of filtered sequences (denoted as N), with the minimum cluster ratio denoted as Rmin and the core cluster ratio denoted as Rcore.

[0073] In this embodiment, Rmin and Rcore are set to 0.015 and 0.03, respectively. These values ​​are non-linearly reduced when N crosses a preset data volume node (such as thousands or tens of thousands) to control the lower limit of the absolute number of clusters, thereby improving the noise resistance and sensitivity of the adaptive balancing algorithm at different sequencing depths.

[0074] (2) Calculate the dynamic threshold according to the formula: N×R, round down, and combine the hard lower limit Min_limit (3 in this embodiment) and the global upper limit Max_limit (200 in this embodiment) for boundary constraints to ensure that the final calculated absolute sequence number is within a reasonable biologically valid range. The formula is: min(Max_limit, max(Min_limit, ⌊N×R⌋)), where ⌊N×R⌋ represents the rounded value of N×R, and R is Rmin or Rcore.

[0075] The core cluster size is a threshold for the number of sequences in the nucleation stage. Because the nucleation stage uses a very high sequence similarity standard, the risk of false positives is low, so a lower absolute number threshold is applied to maximize the capture of potential variant seeds. The minimum cluster size, on the other hand, is the validation threshold for the final consensus sequence output, requiring sufficient sequence redundancy to accumulate within the cluster to ensure statistical support for correcting sequencing errors. This distinction is made to achieve tiered quality control, thereby filtering out background noise that cannot form a stable consensus.

[0076] In this embodiment, the total number of filtered sequences N is 340, ⌊N×Rmin⌋ is 5, ⌊N×Rcore⌋ is 10, and the minimum cluster size is dynamically calculated to be 5, and the core cluster size is 10.

[0077] 2. Progressive clustering: Clustering is performed using a multi-stage thresholding approach. The initial threshold is relatively high (e.g., the 80th to 85th percentile of similarity scores), and the threshold is gradually decreased (e.g., step count = 3) until a final threshold is reached (e.g., the 60th to 70th percentile). This gradual approach helps to preferentially form high-density core clusters.

[0078] In this embodiment, the initial percentile is set to 80, the final percentile to 65, and the number of steps is 3.

[0079] 3. Cluster formation and allocation At each stage, core points and boundary points are identified for each cluster. After the stage ends, unassigned sequences are attempted to be assigned to existing clusters or to form new clusters based on a low similarity threshold (set to 0.3 in this embodiment).

[0080] 4. Cluster size limit To prevent clusters from becoming too large and causing assembly difficulties, the maximum number of sequences in each cluster is limited (e.g., 200 to 600, set to 400 in this embodiment). During the clustering process, once the number of sequences in a cluster reaches the upper limit, the cluster is marked as locked. Thereafter, sequences that should have been assigned to that cluster are retained in the unassigned pool and diverted to candidate clusters with the second highest similarity, or they automatically aggregate into new clusters in subsequent stages. At the same time, clusters that are too small are filtered out (e.g., minimum size = 3, set to 3 in this embodiment).

[0081] Two clustering methods were used: MDS clustering with a stress value of 188.296 and UMAP clustering with a neighbor count of 15 and a minimum distance of 0.1.

[0082] In this embodiment, both clustering processes resulted in three clusters (cluster 1 to cluster 3), as follows: Figure 3 As shown, the sequences within each cluster are saved as independent FASTA / FASTQ files.

[0083] It is worth noting that, due to the high homology of amplicon variants, the overlap of cluster boundaries in the dimensionality-reduced space after the addition of random noise from third-generation sequencing is a normal and expected phenomenon. This will not hinder subsequent assembly, as clustering is only a preliminary step in the "divide and conquer" strategy to reduce computational complexity. The subsequent Optimal Overlap Graph (BOG) algorithm will automatically isolate and remove erroneously mixed edge noise sequences in the graph structure through precise local sequence alignment and reciprocal edge verification.

[0084] S3, perform intra-cluster assembly of sequences based on the optimal overlap graph. The three clusters generated in step S2 are assembled independently using the optimal overlap graph algorithm to generate a consensus sequence representing the cluster.

[0085] I. Overlap Detection The minimap2 tool (with the preset parameter 'ava-ont' in this embodiment, suitable for long read length alignments) is used to detect the overlap between all sequence pairs within the current cluster. The overlapping records are output in PAF format. This step sets the minimum overlap length (100 bp in this embodiment) and the minimum overlap consistency (0.8 in this embodiment).

[0086] II. Construction of Optimal Overlap Graph 1. Preprocessing and Overlap Classification Traverse all alignment records. If there are two or more independent alignment records between any pair of sequences, it is determined that there is a complex repetition region or alignment algorithm artifact between the pair of sequences. These records are marked as "suspicious" and all are removed to prevent chimeric connections.

[0087] After removing suspicious alignments, the proportion of aligned regions to the total sequence length (i.e., alignment coverage) and the proportion of unaligned regions to the total sequence length are calculated. To overcome the inherent high error rate of end sequencing and adapter pruning residue issues in third-generation long-read sequencing platforms, this step introduces a tolerance threshold (e.g., the 98% coverage threshold and 2% boundary tolerance ratio set in this embodiment), and overlaps are classified according to the following criteria: (1) Containment: If the alignment coverage of the sequence is ≥98%, that is, a maximum of 2% end sequencing error is allowed, the sequence is still considered to be completely contained.

[0088] (2) Dovetail: If the inclusion criteria are not met, but the length of the alignment start or end coordinates from the physical endpoint of the sequence is ≤2% (i.e. the alignment has been extended to the biological end within the tolerance range), it is determined to be dovetail.

[0089] (3) Elimination: If none of the conditions are met, the internal overlapping parts are considered to have no biological splicing significance and are directly downgraded and eliminated.

[0090] 2. Create a node Create a graph node for each sequence in the cluster, record the sequence name and length, and initialize the candidate edge set and containment state of the node.

[0091] 3. Create edge and endpoint assignments Extract all alignment records marked as having overlapping ends and create directed edges. Based on the orientation of the alignment, confirm the physical binding positions of the 5' and 3' ends of each edge on the source and target nodes.

[0092] 4. Includes relational clustering and container node extraction Construct an undirected connected graph for edges marked as containing relationships, divide it into multiple containing clusters, and sort them by a comprehensive priority based on three indicators: "sequence length", "graph connectivity" and "number of times it is contained". Select the sequence with the highest priority as the representative node of the container, and mark the remaining nodes as nodes inside the container and hide them.

[0093] 5. Optimal Edge Scoring and Allocation Traverse all unmasked overlapping edges in the graph. For each graph node, extract all candidate edges at its 5' and 3' ends. The algorithm uses a dual sorting rule to determine the best edges: first, sort them in descending order based on overlap length, and then sort them in ascending order based on alignment error rate. Bind the highest-scoring incoming and outgoing edges to the two ends of the node, respectively.

[0094] 6. Reciprocal edge verification Examine all optimal edges. Only when an edge is both the best outgoing edge of the source node and the best incoming edge of the target node is it identified as a mutually beneficial edge and applied to subsequent path solving.

[0095] III. Pathfinding 1. Identify the starting point of the path Prioritize strict starting points (nodes without 5-inch entry edges but with 3-inch exit edges). If no strict starting point exists, use lenient starting points (active nodes with 3-inch edges that are neither contained nor chimeras).

[0096] 2. Path expansion Starting from the origin, extend the path forward along the reciprocal 3' directional edge. During the extension, skip nodes that have been marked as containing or nested, and avoid reusing nodes (to prevent loops).

[0097] 3. Parallel solution Path finding is performed in parallel using multiple threads (8 threads in this example). When the number of starting points exceeds the number of available threads, tasks are allocated using a batch processing method.

[0098] IV. Unitig Construction 1. Path-based Unitig Construct a path-type unitig from the solved path, merging the sequences represented by all nodes on the path according to their overlap. Specifically, the construction steps are as follows: (1) Path traversal and direction correction: Visit the nodes sequentially along the solved valid path. Extract the edge attributes between adjacent nodes. If the edge is marked as reverse, perform reverse complementary processing on the target sequence first.

[0099] (2) CIGAR parsing and coordinate positioning: Extract the CIGAR string from the overlapping records of the two sequences, and combine it with the matching information to locate the actual end coordinates of the overlapping area on the target sequence.

[0100] (3) Redundancy removal and precise splicing: Based on the calculated end coordinates, the redundant parts of the target sequence that overlap with the previous sequence are directly removed, and only the remaining non-overlapping extensions of the target sequence are hard spliced ​​to the end of the current sequence. This process is repeated to complete the sequence fusion of the entire path.

[0101] 2. Processing of remaining nodes After processing the remaining nodes in the pathfinding algorithm, identify the top-level container node, contained nodes, and independent sequences, and construct container-type and independent unitigs respectively. The specific steps are as follows: (1) Node classification screening: The algorithm first traverses the remaining nodes, filters out "contained nodes" (as depth support for folding), and only retains "top-level container nodes" that contain other nodes and "independent nodes" that have no connection.

[0102] (2) Container-type Unitig construction: The complete sequence of the top-level container node is extracted and output separately to retain the longest representative variant with length differences in the sample, and avoid it being discarded as a redundant branch in the regular BOG. ​​The best container node is selected by evaluating the node length, number of connections and number of times it is contained. (3) Independent Unitig construction: Independent sequences are output directly and separately. This step is a fallback to prevent the omission of amplified fragments that are outside the main graph, have extremely low abundance, but are independent.

[0103] 3. Filtering Unitigs with a length less than the minimum threshold (set to 200bp in this embodiment) are filtered out.

[0104] V. Consensus Sequence Construction For each constructed unitig, the original sequencing reads within the cluster are used to construct and refine consensus sequences.

[0105] Two polishing methods are available: Method 1: Use the Racon tool for polishing (e.g., perform 2 rounds by default).

[0106] Method 2: Use step S5 for polishing (based on depth and comparison quality). The difference is that when polishing the consensus sequence within a group, only the reads within that group are used.

[0107] In this embodiment, majority voting is used for polishing, and the number of iterations is set to 1 round.

[0108] In this embodiment, the three clusters each generated their own intra-group consensus sequence files (a total of 13).

[0109] S4, perform inter-group assembly of intra-group consensus sequences based on the optimal overlap graph. All intra-group consensus sequences from different clusters are reassembled to merge those consensus sequences that were separated during the clustering phase but actually have overlapping relationships (e.g., chimeras or highly similar regions).

[0110] 1. Merge all intra-cluster consensus sequences generated in step S3, and add the cluster identifier of its source to each sequence.

[0111] 2. The merged consensus sequence is assembled using the same optimal overlap graph algorithm as in step S3, but with stricter overlap parameters than intra-group assembly. In this embodiment, the minimum overlap length is set to 300bp and the minimum overlap consistency is set to 0.9 to ensure the accuracy of inter-group merging.

[0112] 3. Rollback mechanism In this embodiment, if four consensus sequences are successfully merged, then there will be a final remaining 10 consensus sequences.

[0113] However, in this embodiment, if the inter-group assembly fails (no overlapping edge is found, no path is solved, or a consensus sequence cannot be constructed), the process will fall back to using the intra-group consensus sequence from step S3 as the final result.

[0114] 4. Filtering The final consensus sequence is filtered based on a set expected length (2442bp) and a ratio threshold (e.g., 0.5), and the inter-group assembly result is output. Segments shorter than this ratio threshold typically represent the following invalid cases: (1) Assembly interruption: Due to insufficient local sequencing depth or complex structure, the path solution is broken prematurely, and only a truncated semi-finished sequence is generated.

[0115] (2) Non-target noise: primer dimers and non-specific short fragment amplification products mixed in the sample.

[0116] S5, prunes and refines the consensus sequence based on the comparison depth. Based on the alignment depth information of the original read segments, the consensus sequence is pruned and refined to remove low-quality, low-depth regions and correct bases.

[0117] 1. Comparison Minimap2 is used to align the 340 raw reads filtered in step S1 to the 13 consensus sequences obtained in step S4, generating an alignment file in BAM format.

[0118] 2. Depth computing Traverse the BAM file and calculate the coverage depth for each base position.

[0119] 3. Deep pruning: Trimming at both ends: Starting from both ends of the sequence, examine inwards and remove regions with a depth lower than "average depth × depth ratio threshold" (e.g., the threshold is 0.1 to 0.2, 0.15 in this example).

[0120] Mid-stage pruning: Processes regions in the middle (not at either end) of the sequence whose depth is below a threshold (e.g., average depth × (0.3~0.5), in this example, average depth × 0.4). Three modes are supported: Marking mode: Marks low-depth regions in the output file without changing the sequence.

[0121] Segmentation mode: The sequence is split into two independent sequences at the midpoint of the low-depth region.

[0122] Deletion mode: Delete the low-depth region from the sequence. This embodiment selects this mode.

[0123] 4. Sequence polishing: Based on the alignment quality score and depth information in the BAM file, base correction is performed at each position.

[0124] The final base at that position is determined from the aligned reads using either majority voting or weighted voting (based on quality score).

[0125] In this embodiment, a weighted voting method is selected to correct the bases.

[0126] 5. Iterative pruning Multiple rounds of trimming can be performed (designed as 1 round in this embodiment), and the sequence after the previous round of trimming and polishing is used as the reference sequence for the next round of comparison, thereby gradually improving the sequence quality.

[0127] S6, perform confidence verification on the consensus sequence. The reliability and confidence of the final assembled and pruned consensus sequence are verified.

[0128] The 340 original read segments from step S1 are then compared again to the final consensus sequence after trimming and refining in step S5, and comparison statistics are calculated.

[0129] Low-quality alignments are filtered based on the minimum alignment quality threshold (set to 0.8 in this embodiment) and the minimum alignment length (set to 100bp in this embodiment).

[0130] Calculate various metrics for each consensus sequence, including: average coverage depth, alignment consistency, alignment quality score, and the number of original reads supporting the consensus sequence.

[0131] Based on the confidence score, the best assembly sequence is selected as the representative sequence of the amplicon.

[0132] The 13 consensus sequences obtained in this embodiment are shown in Table 1. Their average coverage depth is 252.659X, average alignment consistency reaches 99.01%, and the number of supported reads is greater than 100, demonstrating the completeness and high accuracy of the assembly results. The consensus sequence with the highest confidence score (unitig_1, 2435bp in length) is almost identical to the validation sequence from the first-generation sequencing: 2338 bases out of 2340 are the same, achieving an alignment consistency of 99.99%, with no gaps.

[0133] Table 1: Statistical Table of Consensus Sequence Results

[0134] This embodiment solves the problems of low efficiency and assembly errors caused by chimeras in existing technologies when processing long read lengths and high similarity amplicon data by using a complete process of sequence filtering, fast clustering based on Strobemer, optimal overlap graph intra-group assembly, inter-group merging, deep pruning, and confidence verification.

[0135] Example 2: Assembly of high-error-rate sequencing data using the method of Example 1. This embodiment applies the method described in Embodiment 1 to ONT shutdown data with a higher error rate. The implementation method is the same as most of the steps in Embodiment 1, with the main difference being the parameter configuration of some modules.

[0136] In this embodiment, a dataset with an expected length of 454bp is processed. After filtering by a factor of 0.6 to 1.5, 41,067 of the original 52,162 read segments remain for assembly. The main difference between this process and Embodiment 1 lies in the following parameters: S2, Sequence clustering based on Strobemer similarity and clustering. In this embodiment, the strobe length k is set to 18, the minimum window offset is set to 10, and the maximum window offset is set to 30.

[0137] S5, prunes and refines the consensus sequence based on the comparison depth. In this embodiment, the regions at both ends are pruned to remove areas whose depth is less than 0.10 times the average coverage depth of the consensus sequence.

[0138] The results show that these reads formed 24 clusters, totaling 72 intra-group consensus sequences. After inter-group consensus construction, 70 sequences were merged, ultimately forming 4 inter-group consensus sequences. One of these was filtered out because its length was less than 1 / 4 of the expected length (113bp). The remaining 3 sequences had an average coverage depth of 12755.2X and an average alignment consistency of 99.06%, but the distribution of the number of supporting reads was uneven, as shown in Table 2.

[0139] Table 2: Statistical Table of Consensus Sequence Results

[0140] This embodiment demonstrates that by adjusting the parameters of different modules, the clustering and graph construction-based method can flexibly adapt to data from different sequencing platforms and with different amplicon lengths, exhibiting strong applicability.

[0141] Example 3: System for assembling amplicon third-generation sequencing data based on clustering and graph construction This embodiment provides a system capable of implementing the method of this application, such as... Figure 4 As shown, it includes the following modules: The data input module is used to: obtain third-generation sequencing data; The sequence clustering module, connected to the data input module, is used to: obtain sequence similarity based on third-generation sequencing data, and cluster the sequences using a clustering algorithm to obtain different clusters; The intra-group assembly module, connected to the sequence clustering module, is used to assemble sequences within the same cluster using a graph-based construction method to obtain intra-group consensus sequences. The inter-group assembly module, connected to the intra-group assembly module, is used to merge all intra-group consensus sequences from different clusters and assemble them using the same graph construction method to obtain inter-group consensus sequences.

[0142] Example 4: An improved system for assembling amplicon third-generation sequencing data based on clustering and graph construction. like Figure 5 As shown, based on the system of Example 3, a sequence trimming and polishing module is also included, which is connected to the data input module and the inter-group assembly module, to implement the following steps: Sequence alignment and depth calculation: The sequence before clustering is used as a reference sequence and compared with the inter-group consensus sequence, and the depth of each base position is calculated; Sequence trimming: Trim regions with a depth lower than the fourth preset threshold; Sequence polishing: Correcting the bases at each position based on alignment quality fraction, depth, and / or number of bases.

[0143] It should be understood that those skilled in the art can use the various forms of processes shown above to rearrange, add, or delete steps. For example, the steps described in this application can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this application can be achieved, and this is not limited herein.

Claims

1. A method for assembling third-generation sequencing data based on clustering and graph construction, characterized in that, The third-generation sequencing data is amplicon third-generation sequencing data, including the following steps: Sequence clustering: Based on the third-generation sequencing data, sequence similarity is obtained, and clustering algorithms are used to cluster the sequences to obtain different clusters; Intra-group assembly: Sequences within the same cluster are assembled using a graph-based method to obtain unitigs. For each constructed unitig, consensus sequences are constructed using the original sequencing reads within the cluster to obtain intra-group consensus sequences. Inter-group assembly: All intra-group consensus sequences from different clusters are merged, and then assembled using the same graph construction method to obtain inter-group consensus sequences. In the sequence clustering step, The sequence similarity is calculated based on Strobemer: a Randstrobes feature vector is generated for each sequence; an index of Strobemer features is constructed, and candidate sequence pairs sharing Strobemer features are found; bidirectional similarity is calculated only for the selected candidate sequence pairs. A multi-stage clustering approach was adopted, with the similarity threshold gradually decreasing at different stages. At each stage, core and boundary points were identified for each cluster. The initial threshold was relatively high, set at the 80th to 85th percentile of the similarity score. Clustering was performed in three steps, with the threshold decreasing at each step, until the final threshold was set at the 60th to 70th percentile of the similarity score. Within-group and between-group assembly are performed based on the optimal overlap graph. The optimal overlap graph assembly includes the following steps: Overlap phase: By comparing sequences, the overlapping relationships between each pair are found. All found overlapping relationships are recorded to generate an overlap graph. Each node in the graph represents a sequence, and edges represent the overlap between two sequences. Mapping stage: The graph is processed to find an optimal path to arrange these sequences and reconstruct the original genome sequence.

2. The method according to claim 1, characterized in that, Before the sequence clustering step, the following steps are further included: Sequence filtering: Filter out sequences whose length is greater than a first preset threshold and whose length is less than a second preset threshold. The first preset threshold and the second preset threshold are determined based on the range of all sequence lengths in the third-generation sequencing data.

3. The method according to claim 2, characterized in that, The sequence filtering step also includes filtering out sequences whose length is less than a third preset threshold.

4. The method according to any one of claims 1-3, characterized in that, Following the intra-group assembly and / or inter-group assembly steps, the process also includes trimming and polishing the sequences: Sequence alignment and depth calculation: The sequence before clustering is used as a reference sequence and compared with the consensus sequence within or between groups, and the depth of each base position is calculated. Sequence pruning: Prune regions with a depth lower than the fourth preset threshold; Sequence polishing: Correcting the bases at each position based on alignment quality score, depth, and / or number of bases.

5. The method according to claim 4, characterized in that, The trimming and polishing steps involve multiple rounds, with the sequence modified and polished in the previous round serving as the reference sequence for comparison in the next round.

6. A system for assembling third-generation sequencing data based on clustering and graph construction, characterized in that, The third-generation sequencing data is amplicon third-generation sequencing data, including the following modules: The data input module is used to: obtain the third-generation sequencing data; The sequence clustering module, connected to the data input module, is used to: obtain sequence similarity based on the third-generation sequencing data, and cluster the sequences using a clustering algorithm to obtain different clusters; The intra-group assembly module, connected to the sequence clustering module, is used to: assemble sequences within the same cluster using a graph-based construction method to obtain a Unitig; and for each constructed Unitig, use the original sequencing reads within the cluster to construct a consensus sequence to obtain an intra-group consensus sequence. The inter-group assembly module, connected to the intra-group assembly module, is used to: merge all intra-group consensus sequences from different clusters, and assemble them using the same graph construction method to obtain inter-group consensus sequences. In the sequence clustering module, The sequence similarity is calculated based on Strobemer: a Randstrobes feature vector is generated for each sequence; an index of Strobemer features is constructed, and candidate sequence pairs sharing Strobemer features are found; bidirectional similarity is calculated only for the selected candidate pairs. A multi-stage clustering approach was adopted, with the similarity threshold gradually decreasing at different stages. At each stage, core and boundary points were identified for each cluster. The initial threshold was relatively high, set at the 80th to 85th percentile of the similarity score. Clustering was performed in three steps, with the threshold decreasing at each step, until the final threshold was set at the 60th to 70th percentile of the similarity score. In both the intra-group assembly module and the inter-group assembly module, assembly is performed based on the optimal overlap graph. The optimal overlap graph assembly includes the following steps: Overlap phase: By comparing sequences, the overlapping relationships between each pair are found. All found overlapping relationships are recorded to generate an overlap graph. Each node in the graph represents a sequence, and edges represent the overlap between two sequences. Mapping stage: The graph is processed to find an optimal path to arrange these sequences and reconstruct the original genome sequence.

7. The system according to claim 6, characterized in that, It also includes a sequence trimming and polishing module, connected to the data input module, and also connected to the intra-group assembly module and / or inter-group assembly module, for implementing the following steps: Sequence alignment and depth calculation: The sequence before clustering is used as a reference sequence and compared with the consensus sequence within or between groups, and the depth of each base position is calculated. Sequence trimming: Trim regions with a depth lower than the fourth preset threshold; Sequence polishing: Correcting the bases at each position based on alignment quality score, depth, and / or number of bases.

8. An electronic device, characterized in that, include: At least one processor; as well as A memory communicatively connected to the at least one processor; wherein, The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform the method of any one of claims 1-5.

9. A non-transitory computer-readable storage medium storing computer instructions, characterized in that, The computer instructions are used to cause the computer to perform the method according to any one of claims 1-5.

Images