A high-throughput single-cell genome collaborative assembly method and system for strain-level resolution

Abstract

A high-throughput single-cell genome collaborative assembly method and system for strain-level resolution. The method first independently assembles and quality assesses a plurality of single-cell amplified genomes (SAGs); then calculates the similarity between SAGs based on the MinHash algorithm, and performs hierarchical clustering to form a co-assembly grouping; for each grouping, an iterative optimization strategy based on tetranucleotide frequency (TNF) analysis is used to identify and eliminate abnormal SAGs; and the optimized SAG combination is used for co-assembly to generate high-quality co-assembled genomes. The application also integrates an automated interactive report generation module, effectively overcoming the defects of uneven coverage, low integrity and severe fragmentation of single-cell genomes, and can distinguish between closely related genomes that are difficult to distinguish in metagenome assembly from mixed samples, providing a high-quality genomic foundation for subsequent strain-level functional gene analysis, and providing a repeatable, high-throughput automated computing framework for microbiome research.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and microbial genomics technology, specifically relating to an automated collaborative assembly and analysis method, system and application for high-throughput single-cell amplified genome sequencing data and related species heterogeneity. Background Technology

[0002] Single-cell genomics and metagenomics technologies provide important tools for elucidating the composition and function of microbial communities, but existing methods still have inherent limitations in handling low-abundance species, distinguishing strain-level differences, and reconstructing uncultured microbial genomes.

[0003] Metagenomics allows us to characterize the composition and function of the microbiome at the community level; however, this method has several inherent limitations. The genomic information of low-abundance species is often masked by signals from high-abundance species, making it difficult to fully interpret. Furthermore, because sequencing data comes from mixtures of multiple species, strain-level genetic variation is usually averaged, leading to an inability to accurately identify functional differences between different strains of the same species. In addition, over 99% of microbial species cannot be obtained through traditional culture methods; this is known as microbial dark matter. The lack of genomic data for these species severely restricts our in-depth understanding of microbial diversity and functional potential.

[0004] Single-cell genomics offers a solution to these problems. This technology, through microfluidic sorting and genome amplification, can bypass the culturing step and directly obtain genomic data from uncultured microorganisms. The preparation of single-cell amplified genomes (SAGs) relies on multiple displacement amplification (MDA) technology, which can extract genomic information from a single cell. With the development of microfluidic technology, SAG preparation has achieved high throughput, routinely generating hundreds to thousands of SAGs in a single study. However, MDA technology itself has inherent biases: amplification bias leads to extremely uneven genome coverage, the amplification of exogenous DNA may introduce contamination, and typically 30%–70% of the genome region cannot be effectively amplified, resulting in generally low genome integrity and severe fragmentation, greatly limiting the depth of subsequent analysis.

[0005] To overcome the aforementioned problems, the SAG co-assembly strategy has emerged. By merging multiple SAG datasets from phylogenetically similar cells and leveraging the complementarity of different SAGs in coverage, assembly quality can be significantly improved. Early ccSAG workflows demonstrated the effectiveness of this strategy, increasing genome integrity from around 40% to over 85% under specific conditions. However, existing SAG co-assembly schemes still have significant limitations: most processes lack complete automation, and steps such as similarity clustering, quality control, and contamination removal still require manual intervention; clustering parameters are mostly based on empirical settings and lack systematic optimization; contamination identification standards are inconsistent, resulting in insufficient reproducibility and difficulty in efficiently utilizing high-performance computing resources.

[0006] Especially in species-rich contexts, such as the oral microbiome which contains over 700 species, and whose strain-level heterogeneity is closely related to diseases like dental caries and periodontitis, significant differences can exist between different strains of the same species in virulence factors, antibiotic resistance genes, and metabolic capabilities. Traditional metagenomic methods and existing SAG co-assembly procedures are insufficient to effectively distinguish these closely related species-level differences, failing to achieve precise species clustering of functional genes. Therefore, developing a fully automated, standardized, scalable, high-throughput SAG co-assembly method and system for strain-level resolution has significant scientific and applied value for research on the oral microbiome and related disease mechanisms. Summary of the Invention

[0007] This invention aims to address the aforementioned problems in existing technologies. Based on Nextflow's CoSAG-nf workflow, this invention constructs a fully automated SAG co-assembly analysis platform that supports large-scale parallel computing and HPC deployment. This workflow system integrates MinHash similarity clustering, tetranucleotide frequency-driven iterative contamination optimization, and a standardized quality control system. Containerization technology ensures cross-platform computational reproducibility. It provides a high-throughput single-cell genome co-assembly method and system for strain-level resolution.

[0008] A high-throughput single-cell amplified genome co-assembly method for strain-level resolution includes (1) performing independent de novo assembly and assembly quality assessment on multiple single-cell amplified genomes respectively; (2) grouping and clustering the single-cell amplified genomes based on the MinHash algorithm to form co-assembly candidate groups; (3) for each co-assembly candidate group, using an iterative optimization strategy based on tetranucleotide frequency analysis to identify and remove abnormal single-cell amplified genomes within the group; (4) using the optimized single-cell amplified genome combination for co-assembly to generate a co-assembled genome, and performing strain-level heterogeneity analysis and functional gene attribution based on the co-assembled genome.

[0009] Optionally, step (5) is also included: serializing the analysis results into structured data and automatically generating a self-contained interactive hypertext markup language report, which includes quality stratification statistics, quality comparison before and after optimization, cluster visualization and classification annotation results.

[0010] In the step of "De novo assembly of single SAGs and assembly quality assessment of single-cell amplified genomes", the single SAG assembly refers to the independent de novo assembly of the sequencing data of each SAG. The --sc parameter is enabled to activate a dedicated algorithm for single-cell data to handle the extremely uneven coverage distribution caused by multiple permutation amplification. A deBruijn graph is constructed through a multi-k-mer strategy, and the graph is simplified and path is parsed using the coverage information. Finally, the contiguous group sequence is output.

[0011] In one specific implementation, an assembly quality assessment system is used to set quality filtering criteria and assess assembly quality; the assessment is based on a lineage-specific marker gene set and uses a machine learning model to estimate genome integrity and contamination rate; integrity represents the degree of completeness of the assembled genome relative to the expected genome, and contamination rate represents the proportion of sequences from other species.

[0012] In another embodiment, the assembly quality assessment system can also provide strain heterogeneity indicators to quantify multiple copy numbers of the same marker gene to reflect potential contamination or actual genome duplication; assemblies that do not meet the quality filtering criteria are marked as low quality and excluded from subsequent co-assemblies, but are still retained in the final report for reference.

[0013] In another implementation, a classification annotation system is used to align bacterial single-copy marker genes, placing the SAG on a reference phylogenetic tree to provide a classification path and calculate similarity metrics such as average nucleotide identity and alignment coverage.

[0014] In one specific implementation, "grouping and clustering SAGs based on genomic similarity estimated by the MinHash algorithm" refers to identifying groups belonging to the same species or closely related strains from SAGs that have passed quality filtering. The MinHash algorithm is used for fast similarity estimation: the genomic sequence is decomposed using k-mer and a hash function is used to select the subset with the minimum value to compress the genome representation. Jaccard similarity is calculated pairwise for all SAGs that have passed quality filtering to construct a similarity matrix, which is then converted into a distance matrix. Subsequently, agglomerative hierarchical clustering is used for grouping, and the clustering quality is evaluated.

[0015] The agglomerative hierarchical clustering supports three connection strategies: full connection, which defines the distance between clusters as the maximum distance between any two points, suitable for strict homogeneous grouping; average connection, which uses the average distance between all point pairs between clusters to balance the compactness and size of clusters; and single connection, which uses the distance between the two nearest points, suitable for capturing gradual population structures; the default is to use full connection to ensure high similarity of co-assembled groups.

[0016] Clustering quality is evaluated using two metrics: silhouette coefficient, which measures the similarity between a sample and its own cluster relative to its nearest neighbor cluster; a higher value indicates a more reasonable clustering. Cophenotypic correlation coefficient measures the correlation between the original distance matrix and the merging distance of the hierarchical tree, and evaluates the fidelity of the clustering tree to the original data.

[0017] The distance threshold for clustering is dynamically determined based on the species definition criteria.

[0018] In one specific implementation, the iterative co-assembly optimization using tetranucleotide frequency analysis and selective read filtering refers to executing an intelligent contamination removal algorithm based on tetranucleotide frequencies for each co-assembly candidate group. This algorithm first calculates the tetranucleotide frequency spectrum for each SAG, identifies and temporarily excludes outlier SAGs through initial outlier detection, and then enters the iterative optimization phase. Employing a greedy optimization strategy, it calculates differentiation scores, constructs experimental co-assembly configurations, evaluates quality based on a hierarchical decision framework, and sets multiple termination criteria to control the optimization process. Finally, it outputs a co-assembled genome with an optimization trajectory log.

[0019] Specifically, the tetranucleotide frequency analysis includes: calculating the tetranucleotide frequency spectrum for each SAG assembly sequence, merging all contigs, statistically analyzing the 256-dimensional tetranucleotide frequency vector using a sliding window approach, and performing Z-score normalization to eliminate the influence of assembly size.

[0020] Specifically, the initial outlier detection stage includes: calculating the Euclidean distance matrix of the tetranucleotide frequency spectrum of all SAGs within the cluster, performing hierarchical clustering using Ward hierarchical clustering, identifying the main cluster and outliers through dendrogram segmentation, and marking SAGs falling outside the main cluster as preliminary outliers and temporarily excluding them.

[0021] Specifically, the iterative optimization phase includes: calculating the average of the tetranucleotide frequency Euclidean distances between each SAG and all other members in the cluster as a differentiation score, sorting them in descending order of differentiation score, and selecting optimization candidates; constructing experimental co-assembly configurations to temporarily exclude each candidate SAG, with co-assembly using a single-cell mode with the --sc parameter enabled; and obtaining integrity, contamination rate, and strain heterogeneity indices through quality assessment of the co-assembly results of each experimental configuration.

[0022] Specifically, the hierarchical decision-making framework includes: a primary objective, namely, achieving a high-quality standard in terms of completeness and a contamination rate below a threshold; a secondary objective, namely, accepting configurations with a contamination rate below a lenient threshold and selecting the configuration with the highest completeness; and a fallback strategy, namely, selecting the configuration with the highest overall score, wherein the overall score is calculated as the weighted value of completeness minus the contamination rate.

[0023] Specifically, the multi-criteria termination conditions include: successful termination, i.e., stopping optimization when the current configuration has reached the primary objective; degradation detection, i.e., regressing and terminating when the completeness of the current iteration decreases by more than a threshold compared to the previous iteration; convergence judgment, i.e. terminating when the quality of all experimental configurations is no better than the current configuration or the number of remaining SAGs is close to the minimum requirement; and maximum iteration limit, i.e. terminating when the preset maximum number of iterations is reached.

[0024] Specifically, the final output co-assembled genome includes an optimized trajectory log, which records the SAG combinations, quality index changes, and decision-making basis for each iteration.

[0025] In one specific implementation, the output report, which includes assembly statistics and / or classification annotation results, refers to an automatically generated interactive, visual, comprehensive quality report. This report serializes all analysis results into a hierarchical data structure, including SAG metadata, assembly statistics, quality indicators, classification annotations, clustering information, and optimization trajectories.

[0026] Specifically, the report includes a quality stratification module that automatically classifies genomes into three levels—high quality, medium quality, and low quality—according to international standards, and calculates the number and proportion of genomes in each quality level to generate a quality improvement matrix before and after optimization.

[0027] Specifically, the report adopts a self-contained architecture design, with all data and dependent libraries embedded in a single HTML file, without external dependencies, and can be opened offline in any modern browser.

[0028] Specifically, the report supports multi-tab navigation, interactive chart displays (including quality distribution scatter plots, clustering efficiency histograms, pollution reduction trend line charts, similarity heatmaps, etc.), dynamic table filtering and sorting, and batch data export.

[0029] This invention further provides a high-throughput single-cell amplification genome co-assembly method and system for strain-level resolution, comprising the following functional modules:

[0030] The single-cell amplified genome assembly and quality assessment module is used to receive sequencing data of single-cell amplified genomes and perform independent de novo assembly, quality assessment, and quality filtering.

[0031] The co-assembly grouping module based on genome clustering is used to calculate the similarity between single-cell amplified genomes using the MinHash algorithm, perform hierarchical clustering analysis, and output co-assembly candidate groups.

[0032] The iterative co-assembly optimization module based on tetranucleotide frequency analysis is used to perform iterative contamination removal driven by tetranucleotide frequency for each candidate group, generate the optimal single-cell amplified genome combination, and perform co-assembly based on the combination to output the co-assembled genome.

[0033] The interactive visualization and comprehensive quality report generation module is used to serialize analysis results into structured data and automatically generate self-contained interactive hypertext markup language reports.

[0034] The system is developed based on a workflow description language and uses containerization technology to encapsulate all software dependencies to ensure cross-platform computing repeatability, supporting large-scale parallel execution in high-performance computing clusters or cloud computing environments. Attached Figure Description

[0035] Figure 1 A schematic diagram of the overall process for CoSAG-nf single-cell genome co-assembly and optimization.

[0036] Figure 2 This is a statistical graph showing the quality distribution of independently assembled single-cell genomes (SAGs); among them, Figure 2 A is the completeness distribution diagram. Figure 2 B is the pollution distribution map.

[0037] Figure 3 Overview of GTDB-Tk classification annotations for single-cell amplified genomes (SAGs) (total SAGs = 7984). Figure 3 A is a pie chart showing the horizontal composition of phyla. Bacillota accounts for the majority (8.5%), with the remainder distributed across more than 10 phyla such as Pseudomonadota and Fusobacteriota. Figure 3 B is a bar chart showing the status of genus-level annotations. 1538 SAGs (19.3%) received genus-level annotations, while 6446 (80.7%) failed to be annotated at the genus level.

[0038] Figure 4 Distribution of high-quality SAGs connectivity and analysis of the top 10 connectivity. Figure 4 A represents the connection number distribution of 5379 high-quality SAGs (log x-axis), with the dashed line indicating the mean (12.37). Figure 4 B represents the 10 SAGs with the highest number of connections and their corresponding number of connections.

[0039] Figure 5Overview of SAG hierarchical clustering results. Figure 5 A represents the distribution of the number of clusters corresponding to each cluster size, with yellow indicating isolated clusters (size = 1) and cyan indicating effective clusters containing multiple SAGs; Figure 5 B represents the composition of all 1,511 clusters, with Multi-SAG clusters accounting for 92.6% and isolated clusters accounting for 7.4%.

[0040] Figure 6 The graph shows the change in contamination level of highly contaminated clusters before and after iterative optimization based on tetranucleotide frequency (TNF). The graph shows 10 representative clusters with a contamination level higher than 10% before optimization. The red bars represent the initial contamination level after co-assembly, and the green bars represent the contamination level after removing abnormal SAGs through TNF iterative optimization.

[0041] Figure 7 CoSAG genome quality assessment results. Figure 7 A is the CoSAG integrity distribution histogram, with colors distinguishing genome quality levels (green: high quality, blue: medium quality), and dashed lines representing the thresholds for medium quality (≥50%) and high quality (≥90%), respectively; Figure 7 B is the histogram of CoSAG pollution distribution, with the dashed line representing the 5% pollution reference threshold.

[0042] Figure 8 Figure showing the GTDB-Tk system classification annotation results of the CoSAG genome. Figure 8 A is a circular diagram showing the horizontal composition of phyla, covering a total of 11 bacterial phyla; Figure 8 B is a horizontal bar chart composed of genera, showing genera that appear ≥ 2 times. Figure 8 C represents the classification resolution statistics: 114 CoSAGs (69.9%) reached the species level annotation, while 49 (30.1%) could only be resolved to the class level. Figure 8 D represents the average nucleotide identity (ANI) distribution between the species-level annotation CoSAG and the most recent reference genome, with a mean of 96.49%, all exceeding the bacterial species classification threshold (95%, red dashed line). Figure 8 E is a horizontal bar chart of the top 30 most frequently occurring species, showing only CoSAGs with species-level annotations.

[0043] Figure 9 Screenshot of the CoSAG-nf process analysis report interface. After the CoSAG-nf process completes, it automatically generates an interactive HTML quality report, covering three modules: process overview, single SAG quality assessment, and cluster analysis.

[0044] Figure 10A heatmap comparing pairwise ANIs of the CoSAG-nf assembled genome and the hybrid assembled reference genome.

[0045] Figure 11 Functional overlap and origin analysis between the CoSAG genome and the B14 reference genome KEGG KO. Figure 11 A is a Venn diagram of the KEGG KO sets of four genome groups (B14, BM, BP, BB), showing the number and proportion of common and unique KOs among the groups; Figure 11 B is based on 1,120 KOs in B14. The UpSet chart shows the distribution of the number of KOs for each source combination. Figure 11 C is a pie chart showing the classification of functional genes in B14, with only 37.3% being core functional genes shared by all three species clusters.

[0046] Figure 12 Comparison of COG functional categories in the CoSAG genome and B14 and distribution of COG-V defense genes. Figure 12 A is the COG functional category deviation heatmap, with B14 as the baseline (Δ% = CoSAG − B14). The colors from blue to red represent the proportions that are below or above B14, respectively. Figure 12 B shows the presence and deletion of 14 CoSAG-specific COG-V class defense genes in four groups of genomes, grouped by functional modules (cytolysin system, restriction modification system, ABC multidrug efflux pump, other defenses).

[0047] Figure 13 A heatmap showing the differences in KEGG pathways between the CoSAG genome and the B14 reference genome. It displays 56 KEGG pathways that differ across the four genome groups, categorized by biological function into antibiotic synthesis and resistance, bacterial motility and biofilm, host-pathogen interactions, and metabolic functions.

[0048] Beneficial effects

[0049] Compared with the prior art, the present invention has the following significant advantages and beneficial effects:

[0050] 1. Fully Automated and Highly Repeatable: The collaborative assembly process based on Nextflow provided by this invention integrates all steps from raw data to the final report, requiring no manual intervention. Containerization technology encapsulates all software dependencies, ensuring computational repeatability across platforms (high-performance computing clusters, cloud platforms, local servers).

[0051] 2. Significantly Improved Genome Quality: Through multi-SAG collaborative assembly, leveraging the complementarity of different SAGs in coverage, fragmented genomes resulting from single-cell sequencing were integrated into genomes with higher continuity and better integrity. In real-world validation with 7,984 oral SAGs, four high-quality and 159 medium-quality genomes were successfully reconstructed from zero high-quality genomes, with taxonomic annotation covering 11 bacterial phyla and a species-level annotation rate of 69.9%.

[0052] 3. Efficient Contamination Removal and Member Optimization: For the first time, an iterative optimization algorithm based on tetranucleotide frequency analysis is systematically introduced into the SAG co-assembly process. This method can accurately identify and remove abnormal SAGs introduced by differentiated lines, cross-contamination, or annotation errors based on clustering, significantly reducing the contamination rate of high-contamination groups from 11%–31% to 3%–10%, effectively avoiding chimerism of co-assembled products, and forming a complete methodological closed loop with the CheckM2 quality assessment standard.

[0053] 4. Overcoming the resolution bottleneck of metagenomics: This invention can resolve microbiome diversity at the strain (species cluster) level. Experiments have shown that this method can successfully disassemble chimeric genomes generated by traditional metagenomic hybrid assembly into multiple independent species clusters with an average nucleotide identity of less than 95%, revealing "micro-diversity" that metagenomic methods cannot identify, and solving the problem pointed out by Oh et al. in skin metagenomics research that "genera-level compositional similarity masks profound differentiation at the strain level".

[0054] 5. Precise Functional Attribution: Supports precise attribution of functional genes at the strain level. Comparative analysis revealed that up to 62.7% of functional genes (including virulence factors and antibiotic resistance genes) in metagenomic hybrid assembly products originate from mixtures of multiple species clusters. This invention can accurately assign these functional genes back to their original species clusters, avoiding systematic biases in functional analysis. For example, biologically significant virulence genes, such as the cytolysin system, are incorrectly attributed to a single species in hybrid assembly products, while this invention enables precise attribution at the species cluster level.

[0055] 6. User-friendly interactive reports: Automatically generates self-contained HTML interactive reports that integrate complete analysis results and quality control charts. Data can be viewed, filtered, and exported in a browser without installing any software, greatly facilitating data sharing and result interpretation.

[0056] 7. Complete Methodological Validation: This invention has achieved complete methodological validation using real oral sample data. The single SAG independent assembly stage yielded only 10 results with medium-quality, high-quality complete genome deletions, which is highly consistent with the inherent coverage non-uniformity of MDA amplification reported by Rinke et al. The formation of 1,399 effective multi-member clusters (92.6%) validates the feasibility of the MinHash hierarchical clustering strategy. The performance of the TNF iterative optimization module is consistent with the key requirements for single-cell genome quality control, ultimately yielding 163 CoSAG genomes, forming a complete technical closed loop from problem identification to solution. Detailed Implementation

[0057] Example 1: Standardized collection of oral samples, preprocessing and sequencing workflow for single-cell genome sequencing

[0058] In this embodiment, 15 volunteers were recruited to participate in saliva collection. After pooling, single-cell genome sequencing was performed. To obtain the best results in single-cell sequencing, it is necessary to remove host cells as much as possible and obtain viable individual microbial cells. The addition of 25% glycerol can effectively protect the activity of microbial cells. After the samples stored at -80℃ were completely thawed on ice, they were added to an equal volume of resuspension and allowed to stand for 30 minutes. They were then filtered sequentially through organic nylon filter membranes with pore sizes of 50 μm and 10 μm. The filtrate was centrifuged and concentrated to obtain concentrated bacterial cells. The volume of each sample was calculated and adjusted according to the absorbance value, and the bacterial count of the 15 samples collected in the morning was pooled. The purpose of sample pooling is to increase the diversity coverage of oral microorganisms, reduce the influence of individual differences, and thus more comprehensively capture representative strains in the oral microbial community to obtain richer single-cell genome data.

[0059] The mixed samples were sent to Motech Biotechnology Co., Ltd. for quality testing, single-cell library construction, and sequencing. The specific process included: preparation of microbial suspension and detection of viability using Live / Dead staining; lysis of microorganisms within droplets to release genomic DNA; whole-genome amplification using MDA technology; DNA fragmentation within droplets; DNA labeling, capture, and PCR amplification within droplets; purification and enrichment of DNA using magnetic beads after droplet disruption; library construction by adding sequencing adapters and index sequences through PCR amplification; binding of the final library to oligos on flowcells; formation of millions of clusters on flowcells through bridging amplification; DNA synthesis using fluorescent dNTPs, synthesizing only one DNA fragment per cycle; photographing and reading the base sequence of each synthesized DNA fragment; and statistical analysis of fluorescence signals to determine the sequence of the template DNA fragment.

[0060] Example 2: Overall Design of the CoSAG-nf Process

[0061] This embodiment describes in detail the specific implementation of each module in the process.

[0062] CoSAG-nf is an end-to-end computational pipeline developed based on Nextflow v23.04.0 for reconstructing high-quality microbial genomes from single-cell amplified genome (SAG) sequencing data. Figure 1 The workflow adopts a modular design and consists of four core modules: (1) Single SAG assembly and quality assessment, which performs de novo assembly and preliminary quality assessment of single-cell amplified genomes; (2) SAG grouping and clustering, which groups SAGs based on genome similarity estimated by the MinHash algorithm; (3) Iterative co-assembly optimization, which removes contamination using tetranucleotide frequency analysis and selective read filtering; and (4) Result visualization and report generation, which outputs a comprehensive interactive HTML report containing assembly statistics and classification annotation results. The workflow accepts single SAG paired-end sequencing data in FASTQ format as input and outputs high-quality co-assembled genomes (FASTA format), a quality assessment report (including indicators such as integrity and contamination rate), and interactive visualization results (HTML format).

[0063] To ensure cross-platform repeatability, all computational steps are standardized using Singularity containerization technology. Based on a dataflow programming model, this process automatically identifies parallel computation steps and supports large-scale parallel execution in high-performance computing (HPC) clusters (such as Slurm, PBS, and SGE) and cloud computing environments (such as AWS and Google Cloud), ensuring the process's efficiency and scalability.

[0064] 2.1 Module 1: Single-cell amplified genome assembly and quality assessment

[0065] This module first performs independent de novo assembly of sequencing data for each SAG. Assembly uses SPAdes v3.15.5, with the --sc parameter enabled to activate a dedicated algorithm for single-cell data. This algorithm can handle the highly uneven coverage distribution (coverage variation coefficient typically > 100%) caused by Multiple Displacement Amplification (MDA). SPAdes constructs a de Bruijn graph using a multi-k-mer strategy (default parameters k=21, 33, 55), and uses coverage information for graph simplification and path resolution, ultimately outputting contig sequences. Assembly quality assessment uses CheckM2 v1.0.2, a tool that estimates genome completeness and contamination rate based on a set of 5,040 lineage-specific marker genes using a machine learning model. Completeness represents the degree of integrity of the assembled genome relative to the expected genome, while contamination rate represents the proportion of sequences from other species. Furthermore, CheckM provides a strain heterogeneity metric, quantifying multiple copy numbers of the same marker gene, reflecting potential contamination or true genome duplication. Quality filtering criteria were set for completeness and contamination rate. Assemblies that did not meet the criteria were marked as low-quality and excluded from co-assemblies, but were still retained in the final report for user reference. Classification annotation was performed using GTDB-Tk v2.3.0 against the GTDB (Genome Taxonomy Database) R214 database, which contains phylogenetic classification information for over 317,000 genomes. GTDB-Tk placed each SAG on a reference phylogenetic tree by aligning 120 bacterial single-copy marker genes, providing a complete classification path from domain to species, and calculating similarity metrics such as average nucleotide identity and alignment coverage, providing taxonomic prior information for subsequent clustering.

[0066] 2.2 Module 2: Identification of Co-assembly Groups Based on Genome Clustering

[0067] This module aims to identify groups belonging to the same species or closely related lineages from thousands of SAGs to support subsequent co-assembly. Considering the computational complexity of whole-genome alignment (O(n²) time complexity), the MinHash algorithm implemented in Sourmash v4.8.0 is used for fast similarity estimation. MinHash compresses the genome representation by performing k-mer (k=31) decomposition on the genome sequence and using a hash function to select the subset with the minimum value (sketch size=1000), making the sketches of similar genomes similar as well, thus transforming the genome alignment problem into a set similarity problem (Jaccard similarity). Jaccard similarity is calculated pairwise for all SAGs that pass the quality filter, constructing an n×n similarity matrix. The similarity matrix is ​​converted into a distance matrix using (1-similarity), and then agglomerative hierarchical clustering is implemented in SciPy v1.11.0.

[0068] The process supports three linkage strategies: (1) Complete linkage defines the distance between clusters as the maximum distance between any two points, tending to produce compact spherical clusters and maximizing the minimum similarity within clusters, suitable for strict homogeneous grouping; (2) Average linkage uses the average distance between all point pairs between clusters, balancing the compactness and size of clusters; (3) Single linkage uses the distance between the two nearest points, tending to produce chain-like clusters, suitable for capturing gradually changing population structures. The default is to use complete linkage to ensure high similarity of co-assembled groups.

[0069] Clustering quality is evaluated using two metrics: (1) the silhouette coefficient, which measures the similarity between a sample and its own cluster relative to its nearest neighbor cluster, ranging from -1 to 1. A higher value indicates a more reasonable clustering; and (2) the cophenetic correlation coefficient, which measures the correlation between the original distance matrix and the merging distance of the hierarchical tree, assessing the fidelity of the clustering tree to the original data. The distance threshold for clustering is dynamically determined based on the species definition criteria, typically corresponding to an ANI ≥ 95%, which is approximately equivalent to Jaccard ≥ 0.90 in terms of MinHash similarity. Each cluster is considered a co-assembly candidate group, requiring a minimum cluster size ≥ 3 SAGs to ensure sufficient coverage complementarity.

[0070] 2.3 Module 3 Iterative Co-assembly Optimization Based on Tetranucleotide Frequency Analysis

[0071] This module is the core innovation of CoSAG-nf, implementing an intelligent contamination removal algorithm based on tetranucleotide frequency (TNF) analysis. A key challenge in co-assembly is that even after clustering based on genome-wide similarity, clusters may still contain differentiated lines, contaminated SAGs, or misclassified annotations. The introduction of these aberrant SAGs significantly increases the contamination rate of co-assembly. TNF analysis utilizes the following biological principle: the tetranucleotide usage pattern of the genome (the frequency spectrum of 256 possible tetranucleotide combinations) is influenced by DNA replication, repair mechanisms, and selection pressure, exhibiting significant conservation in phylogeny, and detectable differences exist in TNF spectra among different species. The algorithm first calculates the TNF spectrum for the assembled sequence of each SAG: all contigs are merged, a 256-dimensional tetranucleotide frequency vector is statistically analyzed using a sliding window approach, and Z-score normalization is performed to eliminate the influence of assembly size. In the initial outlier detection phase, the Euclidean distance matrix is ​​calculated for the TNF spectra of all SAGs within a cluster, and Ward hierarchical clustering is used for hierarchical clustering. The Ward method minimizes the intra-cluster sum of squares increment, tending to produce clusters of balanced size. Dendrogram segmentation identifies the main cluster (typically containing > 80% of SAGs) and outliers. SAGs falling outside the main cluster are marked as preliminary outliers and temporarily excluded.

[0072] The process then proceeds to an iterative optimization phase, employing a greedy optimization strategy to progressively refine the SAG combinations. For the remaining SAGs in the main cluster, the TNF Euclidean distance between each SAG and all other members within the cluster is calculated, and the average value is taken as the "divergence score" of that SAG, reflecting its degree of deviation from the overall population. SAGs are sorted in descending order of their differentiation scores, and the top K SAGs with the highest differentiation scores (default K=3) are selected as optimization candidates. For each candidate SAG, an experimental co-assembly configuration is constructed: the SAG is temporarily excluded, and the sequencing reads of the remaining SAGs are used for co-assembly. Co-assembly uses the SPAdes --sc mode, which is optimized for metagenomic data and can handle mixed reads from multiple strains, separating contigs from different genomes as much as possible through coverage clustering and graph segmentation strategies. The co-assembly results of each experimental configuration are evaluated using CheckM to obtain three indicators: integrity (C), contamination rate (P), and strain heterogeneity (H).

[0073] The quality assessment adopts a hierarchical decision-making framework: (1) Primary objective: integrity ≥ 90% and contamination rate ≤ 5%, meeting the high-quality genome draft (HQ draft) standard; (2) Secondary objective: if the primary objective cannot be achieved, the configuration with a contamination rate ≤ 10% is accepted, and the configuration with the highest integrity is selected from the configurations that meet this condition, balancing contamination control and genome integrity; (3) Fallback strategy: if none of the above standards are met, the configuration with the highest comprehensive score S = C - 5P is selected, which penalizes contamination (the negative impact of contamination is greater than the positive impact of integrity). In each iteration, all experimental configurations (including the exclusion of different candidate SAGs) are compared with the current configuration (without excluding any SAGs), and the configuration with the best quality is selected as the basis for the next round. If a certain excluded configuration wins, the corresponding SAG is permanently removed, the algorithm enters the next iteration, recalculates the TNF differentiation score of the remaining SAGs, and generates a new candidate list.

[0074] The optimization process is designed with multiple termination criteria to prevent over-optimization: (1) Successful termination: The current configuration achieves the target of integrity ≥ 90% and contamination rate ≤ 5%, which is considered to have obtained the optimal result, and the optimization stops; (2) Degradation detection: If the integrity of the current iteration decreases by more than 10% compared with the previous iteration, it indicates that continuing to remove SAGs will start to lose the real genome content, the optimization is over-optimized, and the process will revert to the previous configuration and terminate; (3) Convergence judgment: If the quality of all experimental configurations is not better than the current configuration, or the number of remaining SAGs is close to the minimum requirement (e.g., ≤ 3), the optimization is considered to have converged, the process terminates, and the current best configuration is output; (4) Maximum iteration limit: A maximum number of iterations (default 10) is set to prevent infinite loops. The final output co-assembled genome is accompanied by an optimization trajectory log, which records the SAG combination, quality index changes, and decision basis for each iteration, ensuring that the process is transparent and traceable.

[0075] 2.4 Module 4 Interactive Visualization and Integrated Quality Reporting

[0076] This module transforms complex multidimensional genomic data into user-friendly visualization output. The reporting system uses JSON (JavaScript Object Notation) as the data exchange format and serializes all analysis results into a hierarchical data structure, including: (1) SAG metadata: sample ID, sequencing statistics (number of reads, total number of bases, average read length); (2) assembly statistics: number of contigs, N50, longest contig length, total assembly size; (3) quality indicators: integrity, contamination rate, strain heterogeneity; (4) classification annotation: GTDB classification path, ANI value, classification confidence; (5) clustering information: cluster ID, cluster size, average similarity within the cluster; (6) optimization trajectory: list of SAGs participating in co-assembly, iteration history, quality improvement statistics.

[0077] The quality stratification module automatically classifies genomes according to the international standard (Minimum Information about a Metagenome-Assembled Genome, MIMAG): (1) High Quality: integrity > 90%, contamination rate < 5%; (2) Medium Quality: integrity ≥ 50%, contamination rate < 10%; (3) Low Quality: integrity < 50% or contamination rate ≥ 10%. The statistics module calculates the number and proportion of genomes in each quality level and generates a quality improvement matrix before and after optimization, quantifying the improvement (median, mean) of co-assembly relative to a single SAG in integrity and contamination rate.

[0078] Interactive visualization is implemented using modern Web technology stack: (1) Chart.js v4.2.0 is used to draw bar charts, pie charts, and scatter plots to display genome quality distribution (two-dimensional scatter plot of integrity vs. contamination rate, color-coded classification information), clustering efficiency (cluster size distribution histogram, silhouette coefficient box plot), and classification diversity (species abundance pie chart at the phylum, class, and order levels); (2) Plotly.js v2.20.0 is used for advanced interactive charts, such as a line chart of contamination reduction trend (showing the contamination rate reduction curve during optimization iteration), a heat map (showing the similarity matrix between SAGs), and a three-dimensional scatter plot (simultaneously showing integrity, contamination rate, and genome size); (3) The DataTables jQuery plugin implements the dynamic sorting, search, and pagination functions of the table, and supports multi-dimensional filtering by quality indicators, classification, cluster ID, etc.; (4) The tabbed interface divides the report into Overview, Individual SAGs, Co-assemblies, and Optimization details. Modules such as Details and Taxonomic Summary allow users to quickly switch between views.

[0079] The report employs a self-contained architecture: all data, JavaScript libraries, and CSS stylesheets are embedded in a single HTML file, with no external dependencies. This ensures the report can be opened offline in any modern browser, unaffected by network conditions or server availability, facilitating long-term archiving and cross-institutional sharing. The report also offers batch export functionality, allowing users to export any filtered tables as CSV files, which can then be imported into Excel or R / Python for further statistical analysis. The entire visualization system follows responsive design principles, automatically adapting to different screen sizes (desktop, tablet, mobile phone) and ensuring accessibility through a color-blind-friendly color scheme (Color Universal Design).

[0080] Example 3 Downstream comparative genomics analysis

[0081] To verify CoSAG-nf's ability to finely resolve highly related lineages in complex oral environments, this embodiment selected the complete Streptococcus mitis genome (number B14, CheckM2 assessment of 100% integrity and 0.02% contamination) annotated with GTDB, obtained through a hybrid assembly of second- and third-generation sequencing data, as the gold standard reference. The genetic similarity between the high-quality single-cell genomes (SAGs) assembled by CoSAG and the reference genome B14 was assessed using FastANI v1.34. This software employs the Local Collinear Block (LCB) algorithm, quantifying genetic differentiation by calculating the average nucleotide identity (ANI) between each pair of genomes. Given that Streptococcus mitis is considered a complex containing multiple closely related species clusters under the GTDB taxonomic framework, this study focuses on analyzing CoSAG assemblies (such as Streptococcus mitis_BP, Streptococcus mitis_BM, and Streptococcus mitis_BB) annotated as different subgroups in GTDB. The aim is to verify whether CoSAG-nf can accurately disassemble closely related lineages that are prone to sequence collapse during conventional metagenomic binding, thereby revealing the "microdiversity" in the oral microecology.

[0082] Subsequently, Prokka v1.14.6 was used for gene prediction and preliminary functional annotation, generating GFF format annotation files and FAA format protein sequences. Deep functional annotation was performed using eggNOG-mapper v2.1.12, mapping predicted protein sequences to the eggNOG database via Diamond alignment to obtain multi-level functional information (such as COG classification, KEGG Orthology (KO), pathway, module, GO term, and EC number). To eliminate the influence of genome size and sequencing completeness, percentage normalization was used to calculate the proportion of genes in each category. After excluding genes with unknown functions (S category) and general functional predictions (R category), the tidyverse and ggplot2 packages in R v4.3.0 were used to visualize and compare the remaining 18 core biological functional categories.

[0083] To further identify CoSAG-specific functional components, unique KO entries existing only in the CoSAG genome were screened based on the complete set of KEGG orthologs (KOs) in the B14 genome. These unique genes were then classified and summarized using the COG classification system, with a focus on analyzing defense mechanism categories (COG-V), including ABC type multiple drug efflux pumps, restriction modification systems, toxin-antitoxin systems, and virulence-related factors. By searching the Preferred_name and Description fields (e.g., “cytolysin”, “hemolysin”, “cyl”) in the eggNOG-mapper annotation results, cytolysin and hemolysin-related genes were systematically identified.

[0084] Finally, the distribution heterogeneity of cytolysin and hemolysin genes in the B14 reference genome and CoSAG subgroups (such as Streptococcus mitis_BP, Streptococcus mitis_BM, and Streptococcus mitis_BB) was quantified using Presence / Absence heatmaps, thereby assessing potential differences in virulence function among different species clusters.

[0085] Example 4: Quality Assessment of Single-Cell Genome Sequencing Data

[0086] This embodiment aims to verify the generally limited quality of genomes assembled independently by single cells by evaluating the quality of single-cell amplified genome sequencing data, thereby illustrating the necessity of co-assembly.

[0087] After single-cell genome sequencing and analysis of a pooled sample from 15 individuals, a total of 51 million raw sequences and 7,984 single-cell amplified genomes (SAGs) were obtained. The basic data quality control results are shown in Table 1. The data in the table reveal the sequencing coverage, effective read ratio, and average number of reads per SAG for each sample.

[0088] Table 1 Basic Data Quality Control Table

[0089] After initial independent assembly of each SAG (7984 in total), their quality was assessed and classified using CheckM2 and GTDB-Tk. The results showed that the quality of the genomes assembled independently in single cells was generally limited, with almost all SAGs having less than 50% integrity, and 4695 SAGs having extremely low integrity levels of 0-10%. Figure 2A). Despite good performance in contamination control (contamination levels of 7594 SAGs were between 0-2%), only 10 SAGs met the medium quality standard. Figure 2 (B) This highlights the widespread issues of uneven coverage and high genome fragmentation in single-cell sequencing, further validating the necessity of co-assembling sequence-related closely related SAGs to improve genome quality.

[0090] In terms of taxonomic annotation, due to limitations in completeness, many SAGs fail to obtain clear classification information in GTDB annotations and are usually classified as Unclassified (NA). Figure 3 A) At the genus level, only 1538 (19.3%) of the 7984 SAGs were clearly annotated, while the remaining 6446 (80.7%) failed to resolve to the genus level (see Figure 3). This is closely related to the inherent uneven genome coverage and low integrity characteristics of single-cell sequencing. However, among the identifiable sequences, core phyla of oral microorganisms such as Bacillota, Pseudomonadota, and Fusobacteriota, as well as several functional groups such as Campylobacterota and Bacteroidota, were still accurately detected, presenting a diversity pattern consistent with the expectations of oral microecology. These phyla occupy an important position in oral microecology, demonstrating the effectiveness of the sampling and sequencing process, and also indicating that CoSAG-nf can still retain meaningful taxonomic signals under low-quality input data conditions, laying the foundation for subsequent improvement of genome integrity and taxonomic resolution through co-assembly.

[0091] Example 5: Hierarchical Clustering Analysis Based on Similarity

[0092] This embodiment aims to perform pairwise similarity calculations and hierarchical clustering analysis on all high-quality single-cell amplified genomes using the MinHash algorithm, verifying that the method can effectively identify groups belonging to the same species or closely related strains, and providing candidate inputs for subsequent co-assembly.

[0093] The CoSAG-nf process further utilizes the Sourmash tool to calculate pairwise similarity for all SAGs and explores the structural relationships between SAGs through hierarchical clustering. The analysis sets a minimum similarity threshold of 0.05 and a minimum number of connections of 1 to ensure that the genomes within each cluster have a certain degree of similarity.

[0094] Of the 7984 SAGs, the mean number of connections was 7.66, the median was 2, and the maximum was 141. 2605 low-quality SAGs had fewer than the threshold (0 connections) and were excluded from subsequent cluster analysis, ultimately retaining 5379 high-quality SAGs, a retention rate of 67.4%. Among the high-quality SAGs, the mean number of connections increased to 12.37, the median was 5, and the overall distribution showed a significant right skewness. Figure 4 Figure 4A shows that most genomes have low similarity, while a few SAGs have high connectivity, with the largest number of connections reaching 141, suggesting that these genomes may represent dominant taxa or highly conserved genomic regions in the community (Figure 4B).

[0095] Cluster analysis yielded 1511 clusters, of which 1399 were valid clusters containing multiple SAGs (92.6%), and 112 were isolated clusters containing a single SAG (7.4%). Figure 5 B). The cluster size distribution shows a clear right skewness, with small clusters of size = 2 predominating (584 clusters). As the cluster size increases, the number of corresponding clusters decreases rapidly, and clusters of size ≥ 10 account for only a minority. Figure 5 A). Among the 1399 valid clusters, the average cluster size was 3.8, and the average intra-group similarity was 0.1740. The largest cluster contained 23 SAGs, with an average intra-group similarity of 0.0228, suggesting a certain degree of sequence diversity within this cluster. Overall, the cluster size was highly concentrated in small clusters, reflecting the inherent heterogeneity of single-cell amplification data; a few larger clusters may have enriched dominant species or groups with highly conserved genomic regions, providing candidate inputs for subsequent co-assembly. Figure 5 A).

[0096] Example 6: Collaborative Assembly Optimization and Quality Improvement

[0097] This embodiment aims to verify the effect of the iterative optimization strategy based on tetranucleotide frequency on the removal of contamination from highly contaminated clusters, and to evaluate the improvement effect of co-assembly on genome integrity and contamination level compared with single-cell independent assembly.

[0098] For clusters with contamination levels still above 10% after co-assembly, CoSAG-nf further employed an iterative optimization module based on tetranucleotide frequency (TNF). This module identifies and removes abnormal members within the cluster (including differentiated strains, contaminated SAGs, and misassigned SAGs) round by round, achieving contamination removal and quality improvement. Results showed that this module had a consistent and significant decontamination effect on highly contaminated clusters: before optimization, the contamination levels of most clusters were in the range of 11%–31%, and after TNF-driven iterative removal, the overall level decreased to 3%–10%, with several clusters further reduced to below 5%. Figure 6 This indicates that TNF iterative optimization can effectively reduce the risk of mixed contamination introduced by co-assembly based on clustering and grouping, and can serve as a key quality improvement step after co-assembly.

[0099] After completing similarity-based grouping and TNF-driven iterative optimization, we performed co-assembly on each cluster, ultimately obtaining a batch of CoSAGs (co-assembled SAG genomes), and assessed their genome quality using CheckM2. The results showed that co-assembly significantly improved genome integrity, with 4 genomes reaching high-quality standards and another 159 reaching medium-quality standards. Figure 7 A). The integrity of CoSAG is mainly concentrated in the 50-60% range, accounting for the highest proportion. At the same time, a considerable number of genomes still reach 60-80% integrity, and some genomes are further improved to a higher level of 80-100%. This indicates that co-assembly effectively alleviates the genome fragmentation problem caused by uneven coverage in single-cell sequencing and achieves complementarity of coverage information between different SAGs.

[0100] Regarding contamination levels, the co-assembled CoSAG genomes maintained a low overall level, with most genomes exhibiting contamination levels ranging from 2% to 8%, and a subset of genomes having contamination levels controlled at 0% to 2%, demonstrating good purity. Figure 7 B). Combined with the aforementioned iterative optimization process based on tetranucleotide frequency (TNF), it can be observed that co-assembly not only improves genome integrity but also effectively avoids the problem of contamination accumulation, indicating that the strategy of identifying and eliminating abnormal SAGs plays a key role in maintaining assembly accuracy.

[0101] Overall, compared with the independent assembly stage of a single SAG, co-assembly significantly improved genome continuity and information integrity, while keeping contamination levels within an acceptable range. A batch of medium-to-high-quality CoSAG genomes were successfully obtained, providing a reliable data foundation for subsequent phylogenetic analysis and functional annotation.

[0102] Example 7: Microbial Taxonomy Annotation

[0103] Systematic annotation of the co-assembled CoSAG genome was performed using GTDB-Tk (bac120 marker gene set, GTDB r220 version). A total of 163 CoSAG genomes were successfully annotated, covering 11 bacterial phyla and 29 genera. Figure 8 A). At the phylum level, Bacillota is absolutely dominant with 75 CoSAGs (46.0%), while at the genus level, Streptococcus is the main species (59), consistent with the known ecological position of the oral streptococcal community. Figure 8 B. Figure 8 C). Next were *Pseudomonadota* (27 species, 16.6%), mainly from the *Haemophilus* genus (20 species of *Haemophilus_D*, 3 species of *Haemophilus*); *Patescibacteria* (20 species, 12.3%) contained microscopic parasitic fungi such as HOT-345 and *Nanosynbacter*, which are of significant ecological importance; *Fusobacteriota* (17 species, 10.4%) all belonged to the *Fusobacterium* genus, with *F. pseudooperiodonticum* (7 species) as the representative species. Additionally, oral microecological characteristic groups such as *Bacillota_A*, *Bacteroidota*, *Campylobacterota*, and *Spirochaetota* were also detected.

[0104] Regarding species resolution, 114 CoSAGs (69.9%) achieved species-level annotation, while 49 CoSAGs (30.1%) could only be resolved to the class level. The latter mostly belong to groups with extremely simplified genomes and lacking phylogenetic reference sequences, such as Patescibacteria. ANI analysis of the annotated CoSAGs showed an average ANI of 96.49% (range 95.01%–98.96%) with the nearest reference genome, all exceeding the bacterial species classification threshold (95%), confirming the taxonomic accuracy of CoSAG assembly. Figure 8 D). Among the 114 CoSAGs with species-level annotations, the most frequently occurring species were Haemophilus_D parainfluenzae (12) and Fusobacterium pseudoperiodonticum (7), with Streptococcus oralis_BC (6) ranking third; multiple Streptococcus species and HOT-345-related uncultured taxa were also detected, reflecting the diversity of species composition in oral samples. Figure 8 E).

[0105] Compared with the single SAG stage: Compared with the single SAG independent assembly stage, which resulted in a large number of unclassified classifications due to insufficient completeness, the classification resolution was significantly improved after CoSAG-nf co-assembly, realizing hierarchical annotation from phylum to genus to species, which fully verified the effectiveness of the CoSAG-nf process in improving the completeness of taxonomic information.

[0106] Example 8: Report Data Visualization

[0107] To facilitate systematic review and presentation of analysis results at each stage of the workflow, the CoSAG-nf workflow automatically generates a structured, interactive HTML report after completing all analysis modules. This report integrates key results from the initial SAG quality control to the final CoSAG genome classification and annotation in a single webpage file. It can be viewed directly in a browser without additional software, facilitating data sharing and results reporting within the research group. Figure 9 ).

[0108] The report employs a multi-tab layout, primarily comprising the following functional modules: The Overview page centrally displays core process metrics in the form of statistical cards, including the total number of SAGs processed, the number of co-assemblies formed, the number of high-quality and medium-quality CoSAGs, and the final number of usable genomes. Using this experimental data as an example, the breakdown is as follows: 5379 SAGs, 1511 co-assemblies, 4 high-quality CoSAGs, 159 medium-quality CoSAGs, and 163 final-quality genomes. The Single SAG Quality page displays the assembly quality distribution of a single SAG, including statistical charts of integrity and contamination, as well as CheckM2 evaluation results. The Clustering Analysis page visualizes the preliminary assembly results between SAGs using interactive bar charts and scatter plots, allowing for individual viewing of each cluster member. Furthermore, the report integrates GTDB-Tk taxonomic annotation results and phylogenetic background information, as well as the contamination trend in each iteration of the co-assembly optimization process based on tetranucleotide frequencies.

[0109] The automated generation mechanism of this integrated report significantly improves the repeatability and interpretability of the process, making cross-sectional comparisons between different samples and the optimization and evaluation of process parameters more intuitive and efficient. It also provides a complete analytical record for long-term archiving of experimental data and subsequent reference.

[0110] Example 9 Comparative genomics analysis

[0111] 9.1 Genome sequence comparison and ANI analysis

[0112] Comparative genomic analysis was performed on three representative single-cell co-assembled genomes (Streptococcus mitis_BM, Streptococcus mitis_BP, and Streptococcus mitis_BB) obtained by the CoSAG-nf pipeline, using the second- and third-generation hybrid assembled genome (B14, 100% complete, GTDB annotation: Streptococcus mitis_BM) as a reference. Basic information of each genome is shown in Table 2.

[0113] Table 2 Summary of basic genomic information

[0114] B14 Second and third generation hybrid 100 0.02 2,208,419 Streptococcus mitis BM (cluster_137) CoSAG-nf 72.3 1.8 1,743,991 Streptococcus mitis_BM BP (cluster_1065) CoSAG-nf 92.7 2.4 2,202,651 Streptococcus mitis_BP BB (cluster_1154) CoSAG-nf 50.8 2.7 1,609,023 Streptococcus mitis_BB

[0115] Results of all-pair and all-pair comparisons based on FastANI ( Figure 10 The results showed that the ANI values ​​of all three CoSAG genomes and the reference genome B14 were below the bacterial species definition threshold (95%): the ANI of Streptococcus mitis_BM and B14 was 93.1%, Streptococcus mitis_BB was 91.7%, and Streptococcus mitis_BP was only 84.0%. The pairwise ANIs among the three CoSAG genomes were also below 95% (BM vs BB: 89.6%, BM vs BP: 83.3%, BP vs BB: 82.7%). The results show that the reference genome B14 is actually a chimera of multiple species clusters within the *Streptococcus mitis* complex. Metagenomic sequencing assembles sequences from different species clusters into a single genome, while the CoSAG-nf workflow successfully decomposes it into three species-independent genomic units (*Streptococcus mitis_BM*, *Streptococcus mitis_BP*, and *Streptococcus mitis_BB*) at single-cell resolution. The ANI (Anti-Nutrition Scale) among each species cluster is below the species-defining threshold, confirming their species independence. This result indicates that metagenomic assembly based on mixed sequencing suffers from mixed assembly problems when handling samples of closely related species, while CoSAG-nf has a significant advantage at the species resolution level.

[0116] 9.2 Functional gene overlap and origin tracing

[0117] Based on the functional annotation results of eggNOG-mapper, a total of 1,903 genes were annotated in the B14 reference genome, of which 1,120 (58.9%) obtained KEGG ortholog (KO) numbers. The number of annotated genes in the three CoSAG genomes were 1,362 for Streptococcus mitis_BM, 1,746 for Streptococcus mitis_BP, and 1,182 for Streptococcus mitis_BB, with KEGG annotation rates between 36% and 37%, similar to that of B14, indicating that the annotation depth of each genome is comparable. Venn diagram analysis of the KO sets of the four genomes ( Figure 11 A) shows that the four genomes share a total of 418 core KOs (33.9%), and the three CoSAG genomes have a total of 1,160 KO unions, of which 1,048 (90.3%) are shared with B14; B14 has 72 unique KOs (5.8%). The three CoSAG genomes carry an additional 112 functional genes that are completely deleted from B14, suggesting that CoSAG-nf has supplementary value in terms of the breadth of functional gene detection.

[0118] Tracing the species origin of B14 annotation KO ( Figure 11 B. Figure 11 C) shows that of the 1,120 KOs in B14, only 418 (37.3%) are core functional genes shared by the three CoSAG species clusters, truly representing the conserved functions of the Streptococcusmitis complex; the remaining 630 (56.3%) are from mixed splices of different species clusters, including 263 (23.5%) from BM + BP, 130 (11.6%) from BP + BB, and 51 (4.6%) from BM + BB, as well as 121 (10.8%), 39 (3.5%), and 26 (2.3%) genes unique to BP, BM, and BB respectively but included in B14; the remaining 72 (6.4%) are unique to B14, possibly undetectable due to insufficient CoSAG genome integrity. These results indicate that B14 is a chimera at the functional level, and 62.7% of the annotated functional genes do not truly belong to any single microbial individual.

[0119] 9.3 Comparison of COG Functional Categories

[0120] Comparative analysis based on COG functional categories ( Figure 12A) shows that the proportions of the four genomes in core categories such as translation (J), amino acid metabolism (E), and cell wall / membrane structure (M) are highly consistent, reflecting the basic physiological functional framework shared by the Streptococcus mitis complex. Notable differences include: coenzyme metabolism (H class) is significantly lower in BM (-2.2%) and BP (-2.1%) than in B14; carbohydrate metabolism (G class) is significantly lower in BB (-2.1%); inorganic ion transport (P class) is higher in BM (+1.1%) and lower in BP (-2.0%); and defense mechanisms (V class) are highest in BB (+1.7%), consistent with the finding that this species cluster carries more unique defense genes.

[0121] 9.4 Specificity analysis of defense function genes

[0122] Presence / absence analysis of CoSAG-specific COG-V class defense function genes ( Figure 12 B) A total of 14 gene entries were completely missing in B14, which can be summarized into four categories: (1) Cytolysin system: cylA (K11050) was only found in BB, suggesting that this species cluster has potential cytolysinic toxicity; (2) Type I restriction modification system: hsdS (K01154) was distributed in BM, BP and BB, McrBC methylcytosine restriction system and adenine-specific methyltransferase were specifically detected in BB and BP respectively, and restriction endonuclease (res) was only found in BB. (3) ABC type multidrug efflux pumps: lmrA (K18888) is found in BP and BB, and msbA_1 (K18887) is only found in BP. Both are efflux pump subtypes with complete B14 deletion; (4) Other defenses: Abi-like proteins, β-lactamase C class, glycosyltransferase family 2 and VanZ-like proteins are all specifically distributed in BP, suggesting that BP species clusters have strong immune escape and resistance potential.

[0123] 9.5 KEGG Pathway Differential Analysis

[0124] Differential analysis of the KEGG pathway ( Figure 13 A total of 56 pathways that differed across the four genome groups were identified, which can be categorized into the following types according to their biological functions.

[0125] Regarding antibiotic synthesis and resistance, differential pathways included penicillin / cephalosporin biosynthesis (map00311), neomycin / kanamycin biosynthesis (map00524), vancomycin-like antibiotic synthesis (map01055), and platinum-based drug resistance (map01524). These pathways showed differential distributions among different CoSAG species clusters, which corroborated the results of the aforementioned COG-V class defense gene analysis showing that BP and BB carry the unique ABC efflux pump subtype, suggesting that different species clusters have functional differentiation in antibiotic synthesis and resistance mechanisms.

[0126] Regarding bacterial motility and biofilm, flagella assembly (map02040), bacterial chemotaxis (map02030), and biofilm formation (map02025) were detected in the CoSAG genome but were missing or incomplete in B14. This is consistent with the results of COG analysis, in which class N (cell motility) was higher in CoSAG than in B14, indicating that gene loss occurred during the mixed assembly of this functional module.

[0127] Regarding host-pathogen interactions, the differential distribution of infection-related pathways such as bacterial invasion of epithelial cells (map05100) and pertussis (map05133) suggests that different species clusters may interact differently with the host. Cancer and immune-related pathways in the map05xxx series (including map05203, map05205, map05206, map05226, map05231, and map05340) are mostly found in B14 or parts of CoSAG. These pathways typically reflect functional gene fragments related to host cell processes in bacterial genome annotation, and their differential distribution further supports the conclusion that B14 mixed assembly introduces multi-species functional genes.

[0128] Regarding metabolic functions, the differences in pathways such as lipid metabolism (map01040, map00590), coenzyme synthesis (map00785, map00860), carbohydrate metabolism (map00040, map00053), and secondary metabolite biosynthesis (map00405, map00523) are consistent with the inter-cluster bias of class H (coenzyme metabolism) and class G (carbohydrate metabolism) in COG analysis.

[0129] This demonstrates that CoSAG-nf, with its single-cell genome resolution, successfully disassembled functional chimeras generated by metagenomic hybrid assembly in the same oral sample, achieving accurate species cluster-level attribution of functional genes and revealing the unique toxicity, defense, and metabolic functional differences among species clusters. This proves the significant advantages of this process over traditional metagenomic assembly in functional genomics analysis.

Claims

1. A method for high-throughput single-cell genome co-assembly and analysis at the strain-level resolution, characterized in that, Includes the following steps: (1) Independent de novo assembly and assembly quality assessment were performed on multiple single-cell amplified genomes; (2) Based on the genome similarity estimated by the MinHash algorithm, the single-cell amplified genome is grouped and clustered to form co-assembly candidate groups; (3) For each co-assembly candidate group, an iterative optimization strategy based on tetranucleotide frequency analysis is adopted to identify and eliminate abnormal single-cell amplified genomes within the group; (4) The optimized single-cell amplified genome was used for co-assembly to generate a co-assembled genome, and differential analysis and taxonomic annotation of the co-assembled genome at the level of closely related genomes were performed. Optionally, step (5) is also included: serializing the analysis results into structured data and automatically generating a self-contained interactive hypertext markup language report, which includes quality stratification statistics, quality comparison before and after optimization, cluster visualization and classification annotation results.

2. According to the method of claim 1, the independent de novo assembly in step (1) adopts a single-cell mode assembly algorithm with the --sc parameter enabled, constructs a de Bruijn graph through a multi-k-mer strategy, and uses coverage information for graph simplification and path resolution; the assembly quality assessment is based on the lineage-specific marker gene set to estimate the integrity and contamination rate of the genome, and assemblies that do not meet the quality filtering standards are excluded from subsequent steps.

3. The method according to claim 1, wherein the grouping and clustering based on genomic similarity in step (2) specifically includes: The MinHash algorithm was used to quickly estimate the similarity of single-cell amplified genomes that passed the quality filter. Jaccard similarity was calculated for each pair of single-cell amplified genomes and a distance matrix was constructed. Then, agglomerative hierarchical clustering was used to group them. The clustering distance threshold was dynamically determined according to the species definition criteria, and each cluster was regarded as a co-assembly candidate group.

4. The method according to claim 3, wherein the agglomerative hierarchical clustering supports three connection strategies: full connection, suitable for strict homogeneous grouping; average connection, balancing the compactness and size of clusters; and single connection, suitable for capturing gradual population structure; full connection is used by default; and the clustering quality is evaluated by silhouette coefficient and cophenotypic correlation coefficient.

5. The method according to claim 1, wherein the iterative optimization strategy based on tetranucleotide frequency analysis in step (3) includes: First, the tetranucleotide frequency spectrum of each single-cell amplified genome is calculated and standardized. Outlier single-cell amplified genomes are identified and temporarily excluded through initial outlier detection. Then, the iterative optimization stage is entered. A greedy optimization strategy is adopted to calculate the differentiation score of each single-cell amplified genome, construct experimental co-assembly configurations, evaluate the quality based on a hierarchical decision framework, and control the optimization process through multi-criteria termination conditions.

6. The method according to claim 5, wherein the initial outlier detection comprises: The Euclidean distance matrix was calculated for the tetranucleotide frequency spectrum of all single-cell amplified genomes within the cluster. Ward hierarchical clustering was used for hierarchical clustering. The main cluster and outliers were identified by dendritic segmentation. Single-cell amplified genomes falling outside the main cluster were marked as preliminary outliers and temporarily excluded.

7. The method according to claim 5, wherein the iterative optimization stage comprises: The average Euclidean distance of the tetranucleotide frequencies of each single-cell amplified genome to all other members in the cluster was calculated as the differentiation score. Optimization candidates were selected after sorting the candidates in descending order of differentiation scores. For each candidate single-cell amplified genome, an experimental co-assembly configuration was constructed to temporarily exclude that single-cell amplified genome. Co-assembly was performed in single-cell mode with the --sc parameter enabled. The co-assembly results of each experimental configuration were evaluated for integrity, contamination rate, and strain heterogeneity.

8. The method according to claim 5, wherein the hierarchical decision framework comprises: The primary objective is to achieve a high-quality standard of integrity while keeping the contamination rate below a threshold; the secondary objective is to accept configurations with contamination rates below a lenient threshold and select the one with the highest integrity. The fallback strategy is to select the configuration with the highest overall score, where the overall score is calculated as the weighted value of completeness minus contamination rate. The multi-criteria termination conditions include: successful termination, degradation detection, convergence judgment, and maximum iteration limit.

9. A method and system for high-throughput single-cell amplification of genome co-assembly with strain-level resolution, implementing the method as described in any one of claims 1-8, characterized in that, Includes the following functional modules: The single-cell amplified genome assembly and quality assessment module is used to receive sequencing data of single-cell amplified genomes and perform independent de novo assembly, quality assessment, and quality filtering. The co-assembly grouping module based on genome clustering is used to calculate the similarity between single-cell amplified genomes using the MinHash algorithm, perform hierarchical clustering analysis, and output co-assembly candidate groups. The iterative co-assembly optimization module based on tetranucleotide frequency analysis is used to perform iterative contamination removal driven by tetranucleotide frequency for each candidate group, generate the optimal single-cell amplified genome combination, and perform co-assembly based on the combination to output the co-assembled genome. The interactive visualization and comprehensive quality report generation module is used to serialize analysis results into structured data and automatically generate self-contained interactive hypertext markup language reports.

10. The system according to claim 9, wherein the system is developed based on a workflow description language, uses containerization technology to encapsulate all software dependencies to ensure cross-platform computing repeatability, and supports large-scale parallel execution in high-performance computing clusters or cloud computing environments.

Images