Microbial single cell strain identification method based on difference degree feature screening and multi-index clustering quality evaluation

Abstract

The application provides a microbial single-cell strain identification method based on a difference degree feature screening and multi-index clustering quality evaluation. The method can automatically and accurately identify strains in high-noise and high-sparse microbial single-cell sequencing data through differential SNP feature screening and multi-index clustering quality evaluation, and realizes automatic determination of the number of strains.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and relates to the analysis of microbial single-cell sequencing data, specifically to a method for identifying microbial single-cell strains. Background Technology

[0002] A strain is a smaller distinguishing unit than a species. Different strains may carry different genes, thus exhibiting different biological characteristics. Differences at the genetic level among different strains of the same bacterium directly affect its pathogenicity, metabolic capacity, and environmental adaptability. Therefore, strain-level identification is of great significance in disease treatment, industrial fermentation, and ecological research.

[0003] Currently, strain detection primarily relies on culture-based multiple sequence analysis. This technique requires selecting strains from individual microbial communities, culturing them separately, and sequencing them. This not only favors a small number of easily culturable microorganisms but also results in high costs. While metagenomic sequencing eliminates the need for culture, it often struggles to resolve genomic differences at the strain level. Microplate-based single-cell sequencing can generate strain-level genomes, but its throughput is limited by cost. Microbe-seq, a microbial single-cell genome sequencing technology (Wenshan Zheng et al., High-throughput, single-microbe genomics with strain resolution, applied to a human gut microbiome, Science 376, eabm1483, 2022), utilizes droplet microfluidics combined with bioinformatics analysis to obtain genomic information from tens of thousands of single-celled microorganisms in complex microbial communities without culture. This significantly improves both throughput and resolution, providing a new data foundation for strain-level analysis.

[0004] However, microbial single-cell sequencing data possess inherent characteristics that pose challenges to strain identification. The genome coverage of a single-amplified genome (SAG) is typically only 10%–50%, and the whole-genome amplification process introduces a large number of false-positive single nucleotide polymorphisms (SNPs). Furthermore, each SAG only covers a portion of the reference genome, resulting in a highly sparse SNP matrix formed by the SAG and SNP sites. Accurately identifying closely related strains from such noisy and sparse data is a significant challenge in current microbial single-cell data analysis. Specifically, in the feature selection stage, existing methods typically use variation detection quality indicators such as the QUAL score to filter SNP sites. However, these indicators measure the reliability of variation detection and cannot effectively distinguish between informative and non-informative SNPs for strain identification. A large number of SNP sites with identical genotypes or low coverage consume analytical resources without providing any identification value. Furthermore, the number of SNP sites typically far exceeds the number of SAG sites, and conventional dimensionality reduction methods such as principal component analysis (PCA) assume that the data follows a multivariate normal distribution, which is not suitable for categorical SNP data encoded by allele status and deletion status. In the cluster evaluation stage, existing methods usually rely on single indicators such as the silhouette score to evaluate cluster quality. However, the silhouette score mainly measures the degree of separation between clusters and is not sensitive to the presence of heterogeneous samples within clusters, nor can it reflect the ambiguity of boundary samples. The determination of key parameters such as the number of strains largely depends on human experience, lacking automated, data-driven optimization mechanisms. Additionally, when the number of SAG sites for the target species is small, some dimensionality reduction visualization methods that rely on manifold construction cannot construct an effective low-dimensional representation due to insufficient sample size, making it difficult to provide intuitive verification of strain separation results.

[0005] Therefore, how to automatically and accurately identify strains in high-noise, high-sparse microbial single-cell sequencing data is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0006] To address the aforementioned problems in the existing technology, this invention provides a method for automatically and accurately identifying strains in high-noise, high-sparse microbial single-cell sequencing data. Through differential SNP feature screening and multi-index clustering quality assessment, the number of strains can be automatically determined without supervision, eliminating the need for manual pre-setting of the number of strains.

[0007] The technical solution adopted in this invention is: a method for identifying single-celled microbial strains, comprising the following steps:

[0008] S1. Obtain SNP information of multiple SAGs of the target species relative to the reference genome of the target species, and construct an SNP matrix; in the SNP matrix, the coding value of each SNP site in each SAG is: 1 represents the reference allele, -1 represents the alternative allele, and 0 represents the deletion;

[0009] S2. Sort each SNP site in the SNP matrix and retain the top-ranked SNP sites to obtain the feature SNP matrix; wherein, when sample group information is not provided, calculate the disparity of each SNP site, and sort the SNP sites in descending order according to the disparity, the formula for calculating the disparity is:

[0010] ;

[0011] Where n1 is the number of times the SNP site is encoded as 1 in the SAG, n -1 n is the number of values ​​encoded as -1. non-zero To encode the number of non-zero values, n total The total number of SAGs;

[0012] When sample group information is provided, linear regression is performed using the sample group variables as predictor variables and the code for each SNP locus as the response variable, and the regression determination coefficient R is used as the result. 2 Sort the SNP sites in descending order;

[0013] S3. Perform multiple correspondence analysis to reduce the dimensionality of the feature SNP matrix, obtain the principal coordinates of each SAG in the low-dimensional space, and perform hierarchical clustering of the SAG based on the principal coordinates to obtain a clustering tree;

[0014] S4. Divide the clustering tree into different numbers of candidate strains k, and calculate the comprehensive clustering quality score clean_score for each division result. The number of candidate strains with the highest clean_score is taken as the optimal number of strains. The clean_score is the product of purity rate, boundary NNC, and silhouette coefficient. Wherein, the purity rate is the proportion of SAGs with nearest neighbor consistency equal to 1 to the total number of SAGs, and the nearest neighbor consistency of SAG i is the proportion of its K nearest neighbors belonging to the same strain cluster as SAG i. The boundary NNC is the average of the nearest neighbor consistency of SAGs with nearest neighbor consistency less than 1. The silhouette coefficient is calculated based on the average distance of each SAG to its own cluster of SAGs and the average distance of each SAG to its nearest dissimilar cluster.

[0015] S5. Divide the clustering tree according to the optimal number of strains to obtain each strain cluster; for each SNP site of each strain cluster, count the number of SAG codes 1 and -1 within the cluster, and determine the consensus genotype of the strain cluster according to the majority principle;

[0016] S6. Calculate the consistency ratio of each SAG relative to the consensus genotype of each strain cluster, and assign each SAG to the strain cluster with the highest consistency ratio to obtain the strain identification result; the consistency ratio is: the proportion of SNP sites that encode the same site and are not missing in both SAG and consensus genotype, among which at least one is not missing.

[0017] The main reason why microbial single-cell sequencing data is difficult to use directly for strain identification is that the data is simultaneously high in noise and high in sparsity: amplification errors cause a large number of SNP sites to become false positive signals, while low coverage results in missing values ​​throughout the SNP matrix. To address this, this invention provides a method for automatically and accurately identifying strains from high-noise, high-sparse microbial single-cell sequencing data. Through differential SNP feature screening and multi-index clustering quality assessment, it achieves unsupervised automatic determination of the strain count without the need for manual pre-setting of the strain number.

[0018] First, in the feature screening stage, this invention does not use the conventional variation quality score filtering. Instead, it designs the screening criteria based on the information content of the locus for strain identification: for a SNP locus to distinguish SAGs from different strains, both types of alleles must be present simultaneously in a relatively balanced ratio, and the locus must be actually detected in a sufficient number of SAGs. The difference index combines allele balance and non-deletion coverage into a single score, which precisely characterizes both aspects. Furthermore, when the experimental design includes sample group information such as library batches and processing conditions, linear regression is performed on each locus code using sample group variables, and the results are calculated based on the coefficient of determination R. 2The ranking process allows the utilization of prior grouping information. Two ranking modes cover two common scenarios with and without prior information, providing highly informative feature subsets for subsequent analysis. Secondly, the screened SNP matrix is ​​still categorical coded data; this invention employs multiple correspondence analysis (MCA) for dimensionality reduction. MCA is designed for categorical variables, treating the three coding states of each site as category levels. Its mathematical premise aligns with the {-1,0,1} coded SNP matrix, avoiding the assumption of continuous distribution in conventional dimensionality reduction methods, thus providing a low-dimensional sample representation for hierarchical clustering. Subsequently, this invention uses a comprehensive clustering quality score (clean_score) to iterate through the candidate strains and automatically select the best. This score multiplies the purity rate, boundary nearest neighbor consistency (NNC), and silhouette score: the purity rate answers whether outliers are mixed within a cluster, the boundary nearest neighbor consistency reflects the ambiguity of the transition zone between clusters, and the silhouette score measures the degree of separation between classes. Each of the three factors represents a type of clustering defect; deficiencies in any dimension will lower the overall score. This allows for the exclusion of seemingly separate clusters that are internally mixed or have ambiguous boundaries, ensuring that the automatic determination of strain numbers is based on multi-dimensional evidence rather than relying on a single indicator. Finally, after initial grouping through clustering, the grouping results are further corrected using genotypic evidence: consensus genotypes are generated for each strain cluster according to the majority principle, and the consistency ratio of each SAG relative to each consensus genotype is calculated. The SAG is then assigned to the strain with the highest consistency ratio. Clustering provides the global structure, while genotypic verification provides sample-by-sample molecular evidence; the synergy of both effectively corrects misclassification of boundary samples. These steps work together to ensure the accuracy and automation of strain identification on high-noise, high-sparse data. Verified with real data, this invention can automatically determine the number of strains in multiple species and completely assign all SAGs, demonstrating good clustering quality scores and inter-cluster separation, without requiring manual setting of key parameters.

[0019] Further, the method also includes step S7: generating a template vector based on the observed allele frequency distribution; randomly flipping a predetermined proportion of sites in the template vector to generate one or more (preferably 8-12) synthetic strain genotypes; using the consensus genotypes of each strain cluster and the synthetic strain genotypes to form a reference set; calculating the coverage-weighted distance of each SAG relative to each genotype in the reference set, constructing a distance spectrum matrix with SAGs as rows and each genotype in the reference set as columns; standardizing and reducing the dimensionality of the distance spectrum matrix to generate an embedding graph for verifying strain isolation; in the coverage-weighted distance, the weight w of each site... iThe following rule applies: when both the SAG and the reference genotype are not deleted at this locus and their coding is consistent, w i =1; when both are missing at this site, w i =c r c r The coverage weighting coefficient is defined as [0,1] and is determined based on the coverage of the SAG; when the two encodings at this site are inconsistent, w i =0.

[0020] In step S7, the synthetic strain genotype provides a reference point for each SAG in addition to the real strain, enabling each SAG to obtain a distinguishing distance characterization relative to each reference genotype; the coverage weighting coefficient distinguishes the case of simultaneous deletion of two sites from the case of true consistency, avoiding the misinterpretation of deletion as consistency. The resulting secondary embedding map can provide feature space-based auxiliary verification for strain isolation when conventional dimensionality reduction visualization separation is inadequate.

[0021] Furthermore, the method also includes step S8: merging the sequencing sequences of each SAG within the same strain cluster and assembling them to obtain the genome sequence of each strain; evaluating the completeness and contamination of the genome sequence of each strain, and calculating the assembly quality score clean_score. assembly :

[0022] .

[0023] Step S8 extends the strain identification results to the strain-level genome sequence, and applies a higher weight penalty to the contamination level in the assembly quality score, which facilitates the screening of high-quality strain genomes with high integrity and low contamination.

[0024] Specifically, step S2 further includes: filtering SNP sites with a difference score lower than a difference score threshold, wherein the difference score threshold is preferably 0.05; and retaining N SNP sites. keep Determine as follows:

[0025] ;

[0026] Where, N total is the total number of SNP sites before sorting, and f is the feature retention parameter. Sites with a difference below the threshold usually have a severe allele imbalance or insufficient coverage, and do not provide strain identification value. Direct filtering can reduce the feature dimensionality. The feature retention parameter supports both proportional and fixed-number methods to adapt to SNP matrices of different sizes.

[0027] Specifically, in step S2, the determination coefficient R 2 The calculation formula is:

[0028] ;

[0029] Among them, y ij The encoding of the i-th SNP site in the j-th SAG, ŷ ij For the predicted value of linear regression, ȳ i The mean value encoded by the i-th SNP site; retaining the mean value encoded by the coefficient of determination R. 2 The top 30% of SNP sites in descending order. R 2 R reflects the degree to which site encoding is explained by the sample group variables. 2 High loci are highly correlated with sample grouping and are more likely to carry strain-level identification information.

[0030] Furthermore, the method also includes: traversing multiple parameter combinations of a preset coverage threshold and feature retention parameters, performing steps S2 to S4 for each parameter combination, and selecting the parameter combination that maximizes the clean_score as the optimal parameter configuration. The coverage threshold is used to filter SAGs with genome coverage below the threshold. Thus, key parameters such as the number of retained SNPs and the coverage threshold do not need to be manually set; the system automatically selects the configuration that optimizes clustering quality based on data characteristics, achieving data-driven automatic parameter selection.

[0031] Specifically, in step S4, the number of candidate strains k ranges from 2 to 10.

[0032] Specifically, in step S4, the calculation formula for NNC is:

[0033] ;

[0034] Where, N K (i) represents the K nearest SAGs to SAG i, l j l i These are the strain cluster tags for SAG j and SAG i, respectively.

[0035] Specifically, in step S4, the formula for calculating the contour coefficient is:

[0036] ;

[0037] Where a(i) is the average distance between SAG i and its fellow SAGs, b(i) is the average distance between SAG i and each of the nearest heterogeneous SAGs, and the average value of s(i) of all SAGs is taken as the contour coefficient.

[0038] Specifically, in step S4, if the clean_score corresponding to the number of all candidate strains is not higher than a preset scoring threshold, then all SAGs are determined to belong to the same strain; the preset scoring threshold is preferably 0.2. Setting the scoring threshold enables the method to identify cases where the target species actually contains only a single strain, avoiding forced splitting.

[0039] Specifically, in step S5, determining the consensus genotype of the strain cluster according to the majority principle includes: for the i-th SNP site of strain cluster j, counting the number n of SAG codes that are 1 within the cluster. + j,i The number n encoded as -1 - j,i When n + j,i With n - j,i When n is not greater than 1, the consensus genotype code for that locus is 0; when n + j,i With n - j,i If neither value is 0 and the ratio of the larger to the smaller value is less than 9, the consensus genotype code for that locus is 0; otherwise, the consensus genotype code for that locus is the code corresponding to the larger value. Loci with insufficient evidence or similar proportions of the two alleles are judged as indeterminate to avoid noise loci entering the consensus genotype and interfering with subsequent retesting.

[0040] Specifically, step S6 further includes: if the consistency ratio of a SAG relative to the consensus genotype of all strain clusters is not greater than a preset consistency threshold, or the consistency ratio relative to the consensus genotype of two or more strain clusters is greater than a preset consistency threshold, then the allocation of the SAG is determined to be uncertain and it is assigned to the unknown strain group; otherwise, each SAG is assigned to the strain cluster with the highest consistency ratio; the allocation determination result takes precedence over the clustering result of step S5; after the allocation is completed, strain clusters with no more than 5 SAGs are merged into the unknown strain group. The preset consistency threshold is preferably 0.95. SAGs that are highly consistent with two or more strains should not be forcibly assigned; assigning them to the unknown strain group can avoid misclassification; strains supported by too few SAGs lack statistical reliability; merging them into the unknown strain group can improve the robustness of the identification results.

[0041] In a second aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the above-described method for identifying microbial single-cell strains.

[0042] The beneficial effects of this invention are as follows: This invention uses a dual-mode sorting and screening of SNP sites based on difference indices or group information regression, focusing on sites that are informative for strain identification and reducing feature dimensionality and noise interference; it employs multiple correspondence analysis consistent with classification coding data for dimensionality reduction, providing a reliable low-dimensional representation for hierarchical clustering; it automatically determines the number of strains from three dimensions—intra-cluster purity, boundary clarity, and inter-cluster separation—through a comprehensive clustering quality score based on the product of purity rate, boundary clarity, and silhouette coefficient, eliminating the need for manual parameter tuning; and it re-verifies the clustering results using the consensus genotype consistency ratio, improving the accuracy of strain allocation. Real data verification shows that this invention can automatically determine the number of strains and completely allocate SAGs in multiple species, with clustering quality scores reaching 0.74~0.88 and silhouette coefficients exceeding 0.96; even when conventional dimensionality reduction visualization separation is poor, the strain reference distance spectrum based on coverage-weighted distance still provides clear and reliable strain separation verification. Attached Figure Description

[0043] Figure 1 This is an overall flowchart of the microbial single-cell strain identification method of the present invention;

[0044] Figure 2 UMAP diagram showing the division of 98 SAG strains of Streptococcus parasanguinis_F in Example 1 of this invention;

[0045] Figure 3 This is the reference distance spectrum (Feature UMAP) of the Streptococcus parasanguinis_F strain in Example 1 of the present invention.

[0046] Figure 4 This is the SNP similarity tree of Streptococcus parasanguinis_F in Embodiment 1 of the present invention;

[0047] Figure 5 UMAP diagrams of 38 SAG strains of Streptococcus salivarius in Example 2 of this invention;

[0048] Figure 6 This is the strain reference distance spectrum (Feature UMAP) of Streptococcus salivarius in Example 2 of the present invention.

[0049] Figure 7 This is the SNP similarity tree of Streptococcus salivarius in Embodiment 2 of the present invention;

[0050] Figure 8UMAP diagrams were created for the 357 SAG strains of Gemmiger qucibialis in Example 3 of this invention.

[0051] Figure 9 This is the reference distance spectrum (Feature UMAP) of the Gemmiger qucibialis strain in Example 3 of the present invention.

[0052] Figure 10 This is the SNP similarity tree of Gemmiger qucibialis in Embodiment 3 of the present invention;

[0053] Figure 11 This is a standard UMAP image of Limosilactobacillus urinaemulieris in Embodiment 4 of the present invention;

[0054] Figure 12 This is the reference distance spectrum (Feature UMAP) of the Limosilactobacillus urinaemulieris strain in Example 4 of the present invention.

[0055] Figure 13 This is the SNP similarity tree of Limosilactobacillus urinaemulieris in Embodiment 4 of the present invention. Detailed Implementation

[0056] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0057] Combination Figure 1 The overall process shown in this invention, starting from the input SNP matrix, involves SNP feature screening, multiple correspondence analysis for dimensionality reduction, hierarchical clustering, automatic determination of the number of strains, consensus genotype determination, and confidence allocation retesting to obtain strain identification results. Further, it can generate a strain reference distance spectrum for isolation verification, and perform strain genome assembly and quality assessment, ultimately outputting the strain assembly results (which may include strain genome assembly sequences, strain allocation lists, strain quality statistics, and visualization charts, etc.). The specific implementation of each step is as follows.

[0058] Step S1: Obtain SNP information for multiple SAGs of the target species relative to the reference genome of the target species, and construct an SNP matrix. The SNP matrix is ​​generated by the variant detection module based on a coordinate-sorted VCF file, with the row order maintaining the genomic coordinate order (first by contig, then by locus position). In the SNP matrix, the encoding value of each SNP locus in each SAG is: 1 indicates the reference allele (REF), -1 indicates the alternative allele (ALT), and 0 indicates deletion or indeterminate (No Call). Optionally, the input also includes sample group information metadata and the Guide VCF reference strain genotype; when providing the Guide VCF, only SNPs matching the reference strain loci are retained, and the reference strain genotype is added as a reference column.

[0059] Step S2: Sort each SNP site in the SNP matrix and retain the top-ranked SNP sites to obtain the feature SNP matrix. When sample group information is not provided, calculate the disparity score for each row (corresponding to one SNP site) in the SNP matrix:

[0060]

[0061] in, Let n1 be the number of times the SNP site is encoded as 1 in the SAG, and n be the indicator function. -1 n is the number of values ​​encoded as -1. non-zero To encode the number of non-zero values, n total This represents the total number of SAGs. The dissimilarity value ranges from [0, 0.5]: To measure the balance between the two alleles at this locus. This measures the proportion of SAGs that are actually detected at this locus. SNPs with a difference score below the difference score threshold (default 0.05) are considered non-informative and are directly filtered out. The number of SNPs N retained after sorting by difference score in descending order is... keep Determined by the feature preservation parameter f:

[0062]

[0063] Where, N total This represents the total number of SNP loci before ranking. When sample group information is provided, a linear regression is performed using sample group variables (such as library batch, processing conditions) as predictor variables and the code of each SNP locus as the response variable, and the result is calculated based on the coefficient of determination R. 2 Sort in descending order and keep the highest R 2 Site:

[0064]

[0065] Among them, yij The encoding of the i-th SNP site in the j-th SAG, ŷ ij For the predicted value of linear regression, ȳ i This represents the mean value encoded by the i-th SNP site. In this mode, the top 30% of SNP sites are retained by default. All the above sortings are stable, and sites with the same score maintain their original genomic coordinate order.

[0066] Step S3: Perform multiple correspondence analysis (MCA) to reduce the dimensionality of the feature SNP matrix. MCA treats the three coding states (-1, 0, 1) of each SNP site as three category-level indicator variables, and obtains the principal coordinates of each SAG in the low-dimensional space through generalized singular value decomposition. In this embodiment, the number of principal components is n. components Let's take 5, and the number of iterations be n. iter Set the value to 10. In addition to the feature SNP matrix (feature pattern, default), the input to MCA can also be the SAG similarity matrix (similarity pattern). In the SAG similarity matrix, the encoding vectors a, b ∈ {-1, 0, 1} of two SAGs. m The similarity is expressed as the consistency ratio:

[0067]

[0068] Where, n same n is the number of sites that are not missing and have the same encoding on both sides. diff The similarity is calculated as the number of sites that are not missing on at least one side but have inconsistent encodings. Sites that are missing on both sides are not included in the calculation to avoid interference from missing sites on similarity. When two SAGs have no comparable non-missing sites, the similarity is set to 0.5, and the similarity of the diagonal lines of the same SAG is 1. Subsequently, based on the principal coordinates after dimensionality reduction by MCA, agglomerative hierarchical clustering is used to construct a clustering tree. The Ward method is preferred as the linking method, but the average method can also be used.

[0069] Step S4: Divide the clustering tree into different numbers of candidate strains k (k ranges from 2 to 10). Calculate the overall clustering quality score clean_score for each partition. Select the candidate strain with the highest clean_score as the optimal number of strains. clean_score is the product of purity rate, boundary NNC, and silhouette score.

[0070]

[0071] Nearest-neighbor consistency (NNC) is defined as follows: for each SAGi, take the K nearest SAGs (excluding itself, K is 10 in this embodiment) in its nearest neighbor space, and count the proportion of those SAGs that belong to the same strain cluster as SAGi.

[0072]

[0073] Where, N K (i) represents the K nearest SAGs to SAG i, l j l i These are the strain cluster labels for SAG j and SAG i, respectively; nearest neighbors are determined by Euclidean distance in the principal coordinate space after MCA dimensionality reduction.

[0074] The purity rate is the proportion of SAGs with an NNC of 1 among all SAGs.

[0075]

[0076] Purity rate measures the internal purity of a cluster: if a SAG of another strain is mixed into a cluster of strains, the mixed-in strain is in the heterogeneous neighborhood, its NNC is less than 1, and the purity rate decreases accordingly.

[0077] Boundary NNC is the average NNC of the boundary SAG where NNC is less than 1.

[0078]

[0079] Let the boundary SAG set This refers to the near-neighbor SAGs that are not entirely identical. The boundary NNC is the average NNC of these boundary SAGs, measuring the degree of fuzziness of the cluster boundaries. When the boundaries between clusters are clear, there are few boundary SAGs and their NNCs are close to 1 (only slightly lower than the internal SAGs); when there are large transition regions between clusters, many SAGs are located at the boundaries and their NNCs are significantly lower than 1, thus decreasing the boundary NNC. It complements the purity rate, which only counts the proportion of SAGs with NNC=1 and ignores the fuzziness of the boundary SAGs, while boundary NNC specifically quantifies this part.

[0080] The silhouette score is calculated based on the average distance a(i) between each SAG and its fellow SAGs, and the average distance b(i) between each SAG and the nearest heterogeneous SAGs:

[0081]

[0082] In this embodiment, the silhouette coefficient is calculated based on the cophenetic distance matrix, supplemented by the Calinski-Harabasz index and the Davies-Bouldin index as reference indicators. These three factors penalize three types of defects: intra-class confounding, blurred boundaries, and inter-class non-separation. A deficiency in any dimension will lower the total score, thus excluding clustering results that appear separated but are internally confounding or have blurred boundaries. If the clean_score corresponding to the number of all candidate strains is not higher than the preset scoring threshold (0.2), then all SAGs are determined to belong to the same strain.

[0083] For example, suppose 10 SAGs are clustered into 2 clusters, K=3. Eight of these SAGs have all three nearest neighbors of the same class (NNC=1), while the two boundary SAGs each have one dissimilar nearest neighbor (NNC=2 / 3≈0.67). In this case, the purity is 0.8, and the boundary NNC is 0.67. If the silhouette coefficient is 0.60, then clean_score = 0.80 × 0.67 × 0.60 ≈ 0.32. If the NNC of the same two boundary SAGs is only 0.33 (the boundaries are more blurred), then the boundary NNC is 0.33, and clean_score ≈ 0.16. The total score decreases significantly, accurately reflecting the deterioration of the boundary quality.

[0084] Step S5: Horizontally divide the clustering tree according to the optimal number of strains (cut_tree) to obtain each strain cluster. For each SNP site in each strain cluster, count the number n of SAG codes that are 1 within the cluster. + j,i The number n encoded as -1 - j,i The consensus genotype g of this bacterial cluster was determined based on the majority principle. j,i :

[0085]

[0086] That is: when the counts of both alleles do not exceed 1, the evidence is insufficient and the locus is judged as uncertain (0); when both are non-zero but the ratio of the maximum to the minimum count is less than 9, the proportions of the two alleles are close and are also judged as uncertain (0); otherwise, the code corresponding to the larger count is taken.

[0087] Step S6: For each SAG, calculate its relative to the consensus genotype g of each strain cluster. j Consistency ratio:

[0088]

[0089] This means that only at sites where at least one of the SAG and consensus genotype is not missing, the proportion of sites that are identical in coding and not missing in either genotype is counted. If the consistency ratio of a SAG with the consensus genotype of two or more strain clusters is greater than the preset consistency threshold of 0.95, or if the consistency ratio of a SAG with the consensus genotype of two or more strain clusters is greater than 0.95, then the allocation of the SAG is determined to be uncertain and it is assigned to the unknown strain group; otherwise, the SAG is assigned to the strain cluster with the highest consistency ratio. The retest results take precedence over the initial classification results of hierarchical clustering, that is, the final strain assignment of the SAG is based on the genotype retest; after the allocation is completed, strain clusters with no more than 5 SAGs are merged into the unknown strain group to avoid supporting a strain with too few SAGs.

[0090] Step S7: Generate a strain reference distance spectrum. Based on the observed allele frequency distribution, a template vector is generated. A predetermined proportion (10%) of sites in the template vector are randomly flipped to generate one or more (10 in this example) synthetic strain genotypes.

[0091]

[0092] A fixed reference set is constructed using the consensus genotypes of each strain cluster and the genotypes of the synthesized strains; for each SAG, its coverage-weighted distance relative to each genotype in the reference set is calculated:

[0093]

[0094] Among them, c r The coverage weighting coefficient, ranging from [0,1], is dynamically calculated based on the genomic coverage of the SAG (i.e., the proportion of non-deleted sites in the SAG). When both the SAG and the reference genotype are deleted at a certain site, a partial consistency weight is assigned based on coverage, thus distinguishing between true consistency and deletion consistency. A distance spectrum matrix is ​​constructed with the SAG as rows and each genotype in the reference set as columns. After Z-score standardization, this matrix is ​​used as a substitute feature space for dimensionality reduction, generating a secondary embedding map (Feature UMAP) for verifying strain isolation. This distance spectrum does not participate in the strain allocation determination in step S6; it is only used to assist in verification.

[0095] Step S8: Merge the sequencing reads of each SAG within the same strain cluster (which can be further split by sample), assemble using SPAdes, and obtain the genome sequence of each strain; assess the completeness and contamination of the genome sequence of each strain using an evaluation tool (CheckM, CheckM2, or BUSCO, with CheckM as the default; both CheckM and CheckM2 assess completeness and contamination based on single-copy marker genes), and calculate the assembly quality score:

[0096]

[0097] Assembly quality scores impose a higher weighting penalty on contamination levels. Further screening of strain genomes can be conducted based on quality thresholds: high-quality genomes have at least 90% integrity and no more than 5% contamination; medium-quality genomes have at least 50% integrity and no more than 10% contamination.

[0098] Example 1: Identification of Streptococcus parasanguinis_F strains

[0099] In a human gut microbiome research sample, the single-cell strain identification method of this invention was applied to 98 SAGs of the species *Streptococcus parasanguinis*_F after species identification and binning. Since sample group information was not provided, SNP feature screening was performed using a difference-based model. Subsequently, multi-condition automatic parameter search was performed: 24 parameter combinations were explored, including the coverage threshold *cov_hold* (values ​​0, 0.1, 0.2, 0.3) and the feature retention parameter *apply_feature* (values ​​1000, 2000, 3000, 0.1, 0.2, 0.3). Steps S2 to S4 were performed for each parameter combination. The optimal *clean_score* for each parameter combination is shown in Table 1.

[0100] Table 1. Optimal clean_score for each parameter combination in Example 1

[0101] Coverage threshold cov_hold Feature preservation parameter apply_feature Optimal number of strains k Optimal clean_score Profile coefficient 0.2 0.3 2 0.8834 0.9945 0.2 0.1 2 0.8758 0.9946 0.2 0.2 2 0.8725 0.9923 0.1 0.2 2 0.8432 0.9827 0.1 0.3 2 0.8424 0.9828 0.1 0.1 2 0.8255 0.9802 0.2 3000 2 0.8223 0.9807 0.2 2000 2 0.8096 0.9783 0 0.3 2 0.7929 0.9586 0 0.2 2 0.7896 0.9558 0 0.1 2 0.7820 0.9646 0.2 1000 2 0.7726 0.9665 0.1 3000 2 0.7697 0.9599 0 3000 2 0.7377 0.9500 0.1 2000 2 0.7228 0.9326 0.1 1000 3 0.7074 0.9439 0 2000 4 0.6915 0.9115 0 1000 4 0.6682 0.9168 0.3 1000 2 0.3782 0.9823 0.3 0.3 2 0.3777 0.9828 0.3 2000 2 0.3726 0.9827 0.3 0.2 2 0.3722 0.9782 0.3 3000 2 0.3702 0.9795 0.3 0.1 2 0.3682 0.9751

[0102] As shown in Table 1, the optimal clean_score varies significantly under different parameter combinations. The optimal number of strains was judged to be 3 or 4 under some parameter combinations, which is inconsistent with the correct results. The automatic parameter selection framework selected the parameter combination that maximized the clean_score (cov_hold=0.2, apply_feature=0.3), correctly identifying the optimal number of strains as 2. These results indicate that retaining key parameters such as the number of SNPs has a significant impact on the clustering results, and multi-condition automatic parameter selection can replace manual parameter tuning and provide correct judgments.

[0103] Table 2 shows the clustering quality evaluation indicators corresponding to the number of candidate strains k=2~10 under the optimal parameter combination.

[0104] Table 2. Cluster quality assessment of the number of candidate strains under the optimal parameter combination in Example 1

[0105] Number of candidate strains k clean_score Profile coefficient Calinski-Harabasz index Davies-Bouldin Index 2 0.8834 0.9945 5507.1 0.1072 3 0.5667 0.9780 3680.8 0.3018 4 0.4502 0.9647 2497.2 0.4609 5 0.4158 0.9478 1863.4 0.4311 6 0.3361 0.9231 1489.9 0.4198 7 0.3195 0.9154 1232.9 0.3671 8 0.2927 0.9297 1051.7 0.3154 9 0.2256 0.8355 2721.2 0.3453 10 0.1199 0.8259 2438.4 0.4140

[0106] As shown in Table 2, when k=2, the clean_score is 0.8834, significantly higher than the number of other candidate strains. Therefore, the optimal number of strains is automatically determined to be 2. It is worth noting that when k=4, the silhouette coefficient (0.9647) is still at a high level, but its clean_score (0.4502) is significantly lower than that when k=2. This indicates that relying solely on a single indicator such as the silhouette coefficient may overestimate the classification of strains that have decent inter-class separation but are internally impure or have blurred boundaries. The clean_score, which integrates purity, boundary clarity, and separation, can more accurately identify the optimal number of strains.

[0107] After dividing the cluster tree according to the optimal number of strains (2) and performing consensus genotyping and confidence allocation retesting, 98 SAGs were completely assigned to 2 strains (containing 80 and 18 SAGs respectively), with no undetermined SAGs. Figure 2 UMAP diagrams were generated for 98 SAG strains, showing that two strains (marked by points 0 and 1, respectively) were clearly separated in the two-dimensional embedding space. Figure 3 The reference distance spectrum for strains (Feature UMAP, after removing undetermined SAGs) shows the distinguishing distance of each SAG relative to each strain, with each strain clustering separately and having a clear boundary; Figure 4 The SNP similarity tree (after removing undetermined SAGs) clearly shows that the two strain branches are distinct. The three figures corroborate each other, indicating that the strain division results are stable and reliable.

[0108] Example 2: Identification of strains of Streptococcus salivarius

[0109] For 38 SAGs of the species *Streptococcus salivarius* in the same sample, the same method as in Example 1 was performed. The optimal number of strains was automatically determined to be 2, and the 38 SAGs were completely assigned to 2 strains (containing 23 and 15 SAGs respectively), with no undetermined SAGs; the clean_score was 0.740, the silhouette coefficient was 0.966, the Calinski-Harabasz index was 584.1, and the Davies-Bouldin index was 0.209. Figure 5 UMAP diagrams were generated for 38 SAG strains, with each strain forming a cluster in a two-dimensional embedding space; Figure 6 The reference distance spectrum for the strains (Feature UMAP, after removing undetermined SAG) shows clear separation between the two strains; Figure 7 The SNP similarity tree (after removing undetermined SAGs) clearly shows the branches of the two strains. In this example, the number of SAGs is relatively small (38), yet the method of this invention can still reliably complete strain identification.

[0110] Example 3: Identification of strains of Gemmiger qucibialis

[0111] For 357 SAGs of the Gemmiger qucibialis species in the same sample, the same method as in Example 1 was performed. The optimal number of strains was automatically determined to be 2, and the 357 SAGs were completely allocated to 2 strains (containing 156 and 201 SAGs respectively), with no undetermined SAGs; the clean_score was 0.757, the silhouette coefficient was 0.979, the Calinski-Harabasz index was 8416.1, and the Davies-Bouldin index was 0.187. Figure 8 UMAP diagrams were created for 357 SAG strains. Figure 9 The reference distance spectrum for the strain (Feature UMAP, after removing undetermined SAGs) Figure 10 The SNP similarity tree (after removing undetermined SAGs) shows that the two strains are clearly separated with well-defined boundaries. This example demonstrates that the method of the present invention is also applicable to larger sample sizes with hundreds of SAGs.

[0112] The results of strain identification for Examples 1 to 3 are summarized in Table 3.

[0113] Table 3. Summary of strain identification results for the three species

[0114] Species (GTDB) Total SAG Automatically determined number of strains Strain grouping (SAG number) clean_score Profile coefficient CH Index DB Index Streptococcus parasanguinis_F 98 2 80 / 18 0.883 0.995 5507.1 0.107 Streptococcus salivarius 38 2 23 / 15 0.740 0.966 584.1 0.209 Gemmiger qucibialis 357 2 156 / 201 0.757 0.979 8416.1 0.187

[0115] As shown in Table 3, the method of the present invention automatically determined the number of strains in both species of *Streptococcus* and species of *Gemmiger*, and further distinguished each species into two strains, with all SAGs being completely assigned; the clean_score was 0.74~0.88, and the silhouette coefficient exceeded 0.96, indicating that the present invention can resolve strain-level differences that cannot be distinguished by species identification. In each case, the two strains were well separated and had clear boundaries in the SNP feature space.

[0116] Example 4: Validation of the reference distance spectrum of Limosilactobacillus urinaemulieris strains

[0117] For 110 SAGs of the species *Limosilactobacillus urinaemulieris*, the same method as in Example 1 was used to automatically classify them into two strains (containing 89 and 21 SAGs, respectively). For this case, a conventional UMAP (direct dimensionality reduction of the MCA result) and a strain reference distance spectrum (Feature UMAP) were generated, and the contour coefficients of the two-dimensional embedding were calculated. The results are shown in Table 4.

[0118] Table 4. Comparison of conventional UMAP and strain reference distance spectrum in Example 4

[0119] Visualization Two-dimensional profile coefficient Separation performance Standard UMAP 0.49 The two strains were intertwined at the boundary, making separation difficult. Strain reference distance spectrum (Feature UMAP) 0.75 The two strains were clearly separated with distinct boundaries.

[0120] Figure 11 This is a standard UMAP strain classification diagram. It can be seen that the two strains are intertwined at the boundary, making it difficult to distinguish the strain boundary visually. Figure 12 The reference distance spectrum of the strains (Feature UMAP, after removing undetermined SAG) shows that the two strains are clearly separated and have distinct boundaries. Figure 13 The SNP similarity tree for this species (after removing undetermined SAGs) clearly shows the two strain branches. Combined with Table 4, the two-dimensional profile coefficient increased from 0.49 to 0.75 after using the strain reference distance spectrum, indicating that when conventional low-dimensional visualization cannot clearly present strain separation, secondary embedding based on coverage-weighted distance can provide a more intuitive and reliable verification of strain separation.

[0121] In summary, the method of this invention can automatically determine the number of strains and accurately identify strains in high-noise, high-sparse microbial single-cell sequencing data. All SAGs are completely assigned, and the clustering quality score and inter-class separation degree are both good. Moreover, the results can be verified by strain partitioning map, strain reference distance spectrum and SNP similarity tree.

[0122] In one specific embodiment, an apparatus for identifying microbial single-cell strains is also provided, including a processor and a memory; the processor and the memory are connected via a communication bus; wherein, the processor is used to call and execute a program stored in the memory; the memory is used to store the program, the program being used at least to execute the microbial single-cell strain identification method.

[0123] In one embodiment, a computer-readable storage medium is also provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the microbial single-cell strain identification method as described above.

[0124] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope of the present invention.

Claims

1. A method for identifying single-celled microbial strains, characterized in that, Includes the following steps: S1. Obtain SNP information of multiple SAGs of the target species relative to the reference genome of the target species, and construct an SNP matrix; in the SNP matrix, the coding value of each SNP site in each SAG is: 1 represents the reference allele, -1 represents the alternative allele, and 0 represents the deletion; S2. Sort each SNP site in the SNP matrix and retain the top-ranked SNP sites to obtain the feature SNP matrix; wherein, when sample group information is not provided, calculate the disparity of each SNP site, and sort the SNP sites in descending order according to the disparity, the formula for calculating the disparity is: ； Where n1 is the number of times the SNP site is encoded as 1 in the SAG, n -1 n is the number of values ​​encoded as -1. non-zero To encode the number of non-zero values, n total The total number of SAGs; When sample group information is provided, linear regression is performed using the sample group variables as predictor variables and the code for each SNP locus as the response variable, and the regression determination coefficient R is used as the result. 2 Sort the SNP sites in descending order; S3. Perform multiple correspondence analysis to reduce the dimensionality of the feature SNP matrix, obtain the principal coordinates of each SAG in the low-dimensional space, and perform hierarchical clustering of the SAG based on the principal coordinates to obtain a clustering tree; S4. Divide the clustering tree into different numbers of candidate strains k, and calculate the comprehensive clustering quality score clean_score for each division result. The number of candidate strains with the highest clean_score is taken as the optimal number of strains. The clean_score is the product of purity rate, boundary NNC, and silhouette coefficient. Wherein, the purity rate is the proportion of SAGs with nearest neighbor consistency equal to 1 to the total number of SAGs, and the nearest neighbor consistency of SAG i is the proportion of its K nearest neighbors belonging to the same strain cluster as SAG i. The boundary NNC is the average of the nearest neighbor consistency of SAGs with nearest neighbor consistency less than 1. The silhouette coefficient is calculated based on the average distance of each SAG to its own cluster of SAGs and the average distance of each SAG to its nearest dissimilar cluster. S5. Divide the clustering tree according to the optimal number of strains to obtain each strain cluster; for each SNP site of each strain cluster, count the number of SAG codes 1 and -1 within the cluster, and determine the consensus genotype of the strain cluster according to the majority principle; S6. Calculate the consistency ratio of each SAG relative to the consensus genotype of each strain cluster, and assign each SAG to the strain cluster with the highest consistency ratio to obtain the strain identification result; the consistency ratio is: the proportion of SNP sites that encode the same site and are not missing in both SAG and consensus genotype, among which at least one is not missing.

2. The method for identifying single-celled microbial strains according to claim 1, characterized in that, It also includes step S7: generating a template vector based on the observed allele frequency distribution, randomly flipping a preset proportion of sites in the template vector to generate one or more synthetic strain genotypes; The consensus genotypes of each strain cluster and the genotypes of the synthesized strains together constitute a reference set; Calculate the coverage-weighted distance of each SAG relative to each genotype in the reference set, and construct a distance spectrum matrix with SAG as rows and each genotype in the reference set as columns; The distance spectrum matrix is ​​standardized and dimensionality reduced to generate an embedding graph for verifying strain isolation; in the coverage-weighted distance, the weight w of each site is... i The following rule applies: when both the SAG and the reference genotype are not deleted at this locus and their coding is consistent, w i =1; when both are missing at this site, w i =c r c r The coverage weighting coefficient is defined as [0,1] and is determined based on the coverage of the SAG; when the two encodings at this site are inconsistent, w i =0.

3. The method for identifying single-celled microbial strains according to claim 1 or 2, characterized in that, The process also includes step S8: merging the sequencing sequences of each SAG within the same strain cluster and assembling them to obtain the genome sequence of each strain; assessing the completeness and contamination of the genome sequence of each strain, and calculating the assembly quality score (clean_score). assembly : 。 4. The method for identifying single-celled microbial strains according to claim 1, characterized in that, Step S2 further includes: filtering SNP sites with a difference score below the difference score threshold; and retaining N SNP sites. keep Determine as follows: ； Where, N total denoted as the total number of SNP sites before sorting, and f is the feature retention parameter.

5. The method for identifying single-celled microbial strains according to claim 1, characterized in that, In step S2, the determination coefficient R 2 The calculation formula is: ； Among them, y ij The encoding of the i-th SNP site in the j-th SAG, ŷ ij For the predicted value of linear regression, ȳ i The mean value encoded by the i-th SNP site.

6. The method for identifying single-celled microbial strains according to claim 1, characterized in that, Also includes: Iterate through multiple parameter combinations of the preset coverage threshold and feature retention parameters, and perform steps S2 to S4 for each parameter combination. Take the parameter combination that maximizes clean_score as the optimal parameter configuration.

7. The method for identifying single-celled microbial strains according to claim 1, characterized in that, In step S4, The number of candidate strains k ranges from 2 to 10; The formula for calculating NNC is: ； Where, N K (i) represents the K nearest SAGs to SAG i, l j l i These are the strain cluster tags for SAG j and SAG i, respectively; The formula for calculating the contour coefficient is as follows: ； Where a(i) is the average distance between SAG i and its fellow SAGs, b(i) is the average distance between SAG i and each of the nearest heterogeneous SAGs, and the average value of s(i) of all SAGs is taken as the contour coefficient. If the clean_score corresponding to the number of all candidate strains is not higher than the preset scoring threshold, then all SAGs are determined to belong to the same strain.

8. The method for identifying single-celled microbial strains according to claim 1, characterized in that, In step S5, determining the consensus genotype of the strain cluster based on the majority principle includes: for the i-th SNP site of strain cluster j, counting the number n of SAG codes that are 1 within the cluster. + j,i The number n encoded as -1 - j,i When n + j,i With n - j,i When n is not greater than 1, the consensus genotype code for that locus is 0; when n + j,i With n - j,i If neither value is 0 and the ratio of the larger value to the smaller value is less than 9, the consensus genotype code for that locus is 0; otherwise, the consensus genotype code for that locus is the code corresponding to the larger value.

9. The method for identifying single-celled microbial strains according to claim 1, characterized in that, Step S6 further includes: if the consistency ratio of a SAG relative to the consensus genotype of all strain clusters is not greater than a preset consistency threshold, or the consistency ratio of a SAG relative to the consensus genotype of two or more strain clusters is greater than a preset consistency threshold, then the allocation of the SAG is determined to be uncertain and it is assigned to the unknown strain group; otherwise, each SAG is assigned to the strain cluster with the highest consistency ratio; the allocation determination result takes precedence over the clustering result of step S5; after the allocation is completed, strain clusters with no more than 5 SAGs are merged into the unknown strain group.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the method for identifying microbial single-cell strains as described in any one of claims 1 to 9.

Images