A robust single-cell annotation method based on multiple reference datasets

Abstract

The application discloses a kind of robustness single cell annotation methods based on multiple reference data sets, comprising the following steps: S1: initial reference spectrum construction;S2: initial annotation;S3: confidence evaluation;S4: the construction of local optimization reference spectrum;S5: secondary calibration annotation based on local optimization reference spectrum;S6: result output.The application carries out secondary calibration by introducing local reference spectrum, effectively alleviates the batch effect and heterogeneity difference between reference data and query data, significantly improves the accuracy, robustness and interpretability of single cell annotation.RoMCA method does not need to train special model, has the universality of cross-platform, cross-sample, cross-species, and can be widely used in various single cell annotation tasks, and shows broad application prospect in tumor microenvironment analysis, immunological library research and clinical sample single cell analysis scenes.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and single-cell transcriptome data analysis technology, specifically relating to a robust single-cell annotation method based on multiple reference datasets. Background Technology

[0002] The explosive development of single-cell RNA sequencing (scRNA-seq) technology has made it possible for researchers to analyze gene expression heterogeneity, identify rare cell populations, track developmental trajectories, and characterize the tumor microenvironment at single-cell resolution. However, accurately, robustly, and reproducibly assigning biologically meaningful type labels to each cell, i.e., cell annotation, has always been one of the core bottlenecks between data production and biological discovery.

[0003] Existing annotation strategies can be broadly categorized into two types: Supervised single-reference mapping methods, represented by Seuratv3 / 4, SciBet, SingleR, and scmap, share the common idea of ​​mapping query data to a single, labeled reference dataset. Labels are then directly transferred using similarity scoring or machine learning models. These methods are simple to implement and computationally efficient, but their performance drops sharply when dealing with cross-platform, cross-sample, and cross-species query data. The reasons for this are as follows: The batch effect was not adequately modeled, resulting in distorted similarity measurements, systematic amplification of noise or labeling errors in the reference data, and the algorithm still forcibly assigning meanings when the query data contains new subtypes not covered in the reference, leading to over-annotation and biological false positives. Unsupervised clustering followed by manual matching: For cell type annotation algorithms that use manual matching after unsupervised clustering, systematic bias may occur when dealing with scenarios where there is significant heterogeneity between the reference and query (such as different sequencing platforms, different species, and different disease states) due to the following reasons: The gene detection depth between the 10xChromium and Smart-seq2 platforms can be up to 10-fold, the human-mouse orthologous genes are only about 1:1, and there are a large number of gene family amplification / loss and drug-resistant subgroups or metabolic reprogramming states in the tumor microenvironment that are not covered by the reference dataset. Forced mapping will mask the true biological signals. Summary of the Invention

[0004] (1) Technical problems to be solved To address the shortcomings of existing technologies, the present invention aims to provide a robust single-cell annotation method based on multiple reference datasets. This method addresses the technical problems in existing single-cell transcriptome sequencing data analysis, such as low accuracy, poor robustness, and excessive computational resource consumption in cell type annotation due to differences in sequencing platforms, experimental batch effects, and biological differences between species.

[0005] (2) Technical solution To address the aforementioned technical problems, this invention provides a robust single-cell annotation method based on a multi-reference dataset, comprising the following steps: S1: Initial reference spectrum construction: Obtain single-cell RNA sequencing reference datasets from two or more independent biological samples, and label them with cell type tags. Divide the reference datasets into M reference subsets (M≥3). For each reference subset, group them according to the labeled cell type tags, calculate the average gene expression value of each group of cells, and obtain the pseudo-bulk expression matrix of cell type × gene, which is denoted as the initial reference spectrum matrix. S2: Initial annotation: Extract the gene expression vector of each cell in the single-cell sequencing dataset to be queried. Based on the M initial reference spectrum matrices obtained in step S1, solve the non-negative least squares model to obtain the M initial cell type annotations for each cell. S3: Confidence assessment: For each cell to be queried, count the cell type that appears most frequently in the M initial annotations and calculate the proportion of the occurrences of this type to M; when the proportion is greater than or equal to a preset threshold, the cell is classified into a high-confidence cell subgroup, otherwise it is classified into a low-confidence cell subgroup. S4: Construction of the local optimized reference spectrum: Extract gene expression data and corresponding cell type labels from high-confidence cell subpopulations, recalculate cell type-specific pseudo-bulk expression profiles, and obtain the local optimized reference spectrum matrix; the local optimized reference spectrum matrix relies only on the high-confidence information of the query data itself, has stronger biological consistency and batch adaptability, and outliers are truncated during the standardization process to ensure the stability and reliability of the matrix; S5: Secondary calibration annotation based on local optimized reference spectrum: Using the local optimized reference spectrum matrix obtained in step S4 as the sole benchmark, non-negative least squares analysis is performed again on all cells to be queried to obtain the cell type annotation results of the second annotation; during the secondary annotation process, a verification label is added to cells with low confidence, and subpopulation clustering, differential expression analysis or comparison with external databases can be performed to identify potential new subpopulations. S6: Output Results: The first and second rounds of annotation results are integrated into the final annotation results, confidence scores, and local optimized reference spectrum matrix, and written to a structured file in CSV or rds format to support efficient storage and subsequent analysis; optionally, the annotation results and confidence distribution can be visualized in a low-dimensional embedding space.

[0006] Furthermore, the splitting method must ensure the independence and diversity of the reference subsets. The reference subsets are split according to biological individuals, sequencing batches, or platform sources, and the number of cells in each subset is ≥500. Furthermore, for subsets with fewer than 500 cells, random oversampling or merging into adjacent subsets is used to ensure statistical robustness. The splitting method further includes dividing according to disease state, tissue location, or time point to increase the biological diversity of the reference subsets.

[0007] Furthermore, the objective function of the nonnegative least squares model is min‖y-Xβ‖ 2 2. The constraint condition β≥0, where y is the gene expression vector of a single query cell, X is the reference spectrum matrix, and β is the cell type contribution coefficient vector. After solving the problem, the cell type corresponding to max (β) is used as the annotation result under this reference subset. Normalize β; The model is solved using a low-level R implementation to support efficient parallel computation, accelerated by pbmclapply, and supports computations >10. 4 Rapid processing of large numbers of cells.

[0008] Furthermore, the preset threshold is 0.7-0.9, more preferably 0.8; when the ratio is exactly equal to the threshold, it is defaulted to be classified into a high-confidence cell subpopulation; the preset threshold can be dynamically adjusted according to the batch effect strength, platform differences or species differences among reference subsets, and the dynamic adjustment of the preset threshold should be based on ensuring the purity and representativeness of the high-confidence subpopulation as the core objective.

[0009] Furthermore, after the secondary annotation in step S5 is completed, it further includes: Individual source tracing analysis is performed on low-confidence cell subpopulations, and their expression profiles can be displayed as similarity heatmaps with the local optimized reference spectrum matrix. High-confidence cell subpopulations can be further subdivided into functional subpopulations by using graph clustering or density clustering methods. The source tracing analysis further includes comparing low-confidence cells with external databases (such as CellMarker and PanglaoDB) to identify potential novel cell types or transitional cells.

[0010] Furthermore, a gene screening step is included before step S1: (1) The COSG differential expression algorithm was used to identify cell type-specific marker genes from the merged reference dataset; (2) Perform fold change analysis on the pseudobulk expression profile and retain genes with FoldChange>1.5 among different cell types; (3) Using the gene list obtained in steps S1 and S2, the reference spectrum matrix and the dataset to be queried are synchronously subsetted to reduce noise and improve computational efficiency; the gene screening step further includes removing low-expression genes or genes with high deletion rates to ensure the quality and robustness of the input data.

[0011] A robust single-cell annotation system based on a multi-reference dataset includes: The initial annotation module is used to execute steps S1 to S3 and output high / low confidence cell subpopulations. The local reference building module is used to execute step S4 and generate a local optimized reference spectrum matrix; The secondary annotation module is used to execute step S5 and output the final annotation result; The results visualization module is used to display annotation results and confidence distribution in a low-dimensional space; The logging and auditing module records the parameters, version number, and execution time of each annotation task to ensure that the results are traceable. It also includes a user access control module to support concurrent access by multiple users and data isolation.

[0012] Furthermore, both the initial annotation module and the secondary annotation module are implemented using the same non-negative least squares solver, and the solver supports parallel computation to accelerate processing >10. 4 Processing of cells at the scale.

[0013] Furthermore, the results visualization module supports multiple dimensionality reduction algorithms, including UMAP and t-SNE, and supports interactive selection and export of high-confidence cell subpopulations. The results visualization module further includes a confidence level slider, customizable cell type coloring scheme, and dynamic legend generation to improve user experience and result interpretation efficiency.

[0014] Beneficial effects Compared with the prior art, the beneficial effects of the present invention are as follows: This invention effectively mitigates batch effects and heterogeneity differences between reference and query data by introducing a locally optimized reference spectrum for secondary calibration, significantly improving the accuracy, robustness, and interpretability of single-cell annotation. Furthermore, by constructing a two-stage non-negative least squares iterative analytical framework, this invention achieves adaptive high-precision annotation across datasets without the need for pre-training complex classifiers. The proposed method does not require training a dedicated model and possesses universality across platforms, samples, and species, making it widely applicable to various single-cell annotation tasks. It demonstrates broad application prospects in scenarios such as tumor microenvironment analysis, immune repertoire research, and single-cell analysis of clinical samples. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the single-cell data annotation process of the present invention; Figure 2 This is a performance comparison chart of the RoMCA method and SciBet on the Wilk dataset provided in this embodiment of the invention. Figure 3 A performance comparison chart of the RoMCA method provided in this embodiment of the invention with existing mainstream cell type annotation algorithms across datasets; Figure 4 The ablation experimental results of each module of the RoMCA method provided in the embodiments of the present invention are shown. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0017] It should be noted that the use of terms such as "an embodiment," "embodiment," and "exemplary embodiment" in the specification indicates that the described embodiment may include a specific feature, structure, or characteristic, but not every embodiment necessarily includes that specific feature, structure, or characteristic. Furthermore, when describing a specific feature, structure, or characteristic in conjunction with embodiments, implementing such a feature, structure, or characteristic in conjunction with other embodiments should be within the knowledge of those skilled in the art.

[0018] This invention provides a robust single-cell annotation method based on a multi-reference dataset. This method is executed on a processor via computer program instructions. Through a collaborative design of multi-reference subset construction, parallel model solving, and ensemble confidence assessment, it effectively resists interference from batch effects, sequencing bias, and cell heterogeneity, improving the accuracy and biological reliability of the annotation results. Specifically, it includes the following steps: S1. Initial reference spectrum construction: Obtain single-cell RNA sequencing reference datasets from two or more independent biological samples, and label them with cell type tags. Divide the reference datasets into M reference subsets (M≥3). For each reference subset, group them according to the labeled cell type tags, calculate the average gene expression value of each group of cells, and obtain the pseudo-bulk expression matrix of cell type × gene, which is denoted as the initial reference spectrum matrix.

[0019] The splitting method must ensure the independence and diversity of the reference subsets. The reference subsets are split according to biological individuals, sequencing batches, or platform sources, and the number of cells in each subset is ≥500. Furthermore, for subsets with fewer than 500 cells, random oversampling or merging into adjacent subsets is used to ensure statistical robustness. The splitting method further includes dividing according to disease state, tissue location, or time point to increase the biological diversity of the reference subsets.

[0020] Step S1 aims to establish a reference space with high biological diversity and free from systematic bias for subsequent annotation tasks, providing a solid data foundation for accurate annotation. Its core principle is to balance the coverage and consistency of the reference data, specifically: First, obtain one or more single-cell RNA sequencing (scRNA-seq) reference datasets that have been experimentally validated and have accurate cell type annotations. These reference datasets are preferably from public databases (such as GEO and ArrayExpress) or self-built datasets that have undergone strict quality control to ensure the authority and reliability of cell type annotations. Secondly, in order to fully cover biological variations from different donors, experimental batches, and physiological states, and to accurately simulate multi-reference scenarios in real-world applications, the system logically splits the reference dataset into M independent reference subsets based on key metadata of the sample source (e.g., donor ID, sequencing batch number, experimental processing conditions, tissue source, or species information). M is set to an integer greater than or equal to 3—this lower limit is based on statistical principles, ensuring that subsequent ensemble analysis has a sufficient sample size to reduce random errors, while avoiding the problem of insufficient diversity of reference subsets due to an excessively small M value. Subsequently, for each of the split reference subsets, the system performs pseudo-bulk expression profile calculations. The core purpose of this operation is to address the inherent sparsity problem of scRNA-seq data while preserving expression features at the cell type level. Specifically, the system first preprocesses and filters the single-cell data within the reference subsets, removing low-quality cells (such as cells with an excessively high proportion of mitochondrial genes or abnormal gene detection numbers). Then, it aggregates the gene expression vectors of all qualified single cells belonging to the same cell type label to further improve the reliability of the pseudo-bulk expression profile. Through this aggregation operation, the originally sparse single-cell expression matrix with the dimension of single cell number × feature gene number is transformed into a dense initial reference profile matrix with the dimension of feature gene number × cell type number. Finally, for the M reference subsets, M independent initial reference profile matrices are constructed, each matrix corresponding to the cell type expression features of a set of reference subsets.

[0021] Furthermore, significant systematic technical deviations often exist between different reference subsets due to factors such as sequencing depth, library construction size, cell capture efficiency, and differences in sequencing platforms. These deviations can severely interfere with the accuracy of subsequent similarity calculations, leading to distorted annotation results. Therefore, this invention performs column-direction standardization on the M initial reference spectrum matrices one by one. By uniformly projecting all reference spectra onto a unit feature space hypersphere, not only are the differences in expression magnitude caused by technical deviations eliminated, but the fairness of subsequent similarity calculations based on Euclidean distance, cosine similarity, or projection mapping is also ensured, making the cell type expression features of different reference subsets comparable.

[0022] S2: Initial annotation: Extract the gene expression vector of each cell in the single-cell sequencing dataset to be queried. Based on the M initial reference spectrum matrices obtained in step S1, solve the non-negative least squares model to obtain the M initial cell type annotations for each cell. The objective function of the nonnegative least squares model is min‖y-Xβ‖. 2 2. The constraint β≥0, where y is the gene expression vector of a single query cell, X is the reference spectrum matrix, and β is the cell type contribution coefficient vector. After solving the problem, the cell type corresponding to max(β) is used as the annotation result under this reference subset; further, β is normalized. The model is solved using a low-level R implementation to support efficient parallel computation, accelerated by pbmclapply, and supports computations >10. 4 Rapid processing of large numbers of cells.

[0023] Step S2 utilizes the diverse reference benchmarks constructed in Step S1 to perform a preliminary global exploration of the test data. The goal is to quickly obtain candidate cell type labels for each query cell and provide multi-perspective results support for subsequent confidence assessment. Specifically: First, the single-cell dataset to be queried is read, and the same preprocessing procedure as the reference dataset is performed on it (including low-quality cell filtering, gene expression normalization, and feature gene screening) to ensure that the data format and quality are consistent. Then, for each query cell, its high-dimensional gene expression vector is extracted as the input for subsequent model solving. Based on the M standardized initial reference spectrum matrices constructed in step S1, M non-negative least squares (NNLS) models are independently constructed and solved for each query cell. A parallel computing architecture is used to allocate computing resources to avoid efficiency bottlenecks caused by a large value of M or a large number of query cells. The core mathematical expression of the model is: under the non-negative constraint β≥0, find the optimal non-negative coefficient vector β such that the objective function min‖y-Xβ‖ 2 2. Minimize, where y is the gene expression vector of a single query cell, X is a single normalized initial reference spectrum matrix (dimension is the number of feature genes × number of cell types), and β is the cell type contribution coefficient vector (dimension is the same as the number of cell types, and each element corresponds to the contribution ratio of a cell type to the expression feature of the query cell). Introducing a non-negative constraint of β≥0 perfectly aligns with the biological essence—gene expression level, as an indicator reflecting the degree of gene activation, only differs in whether it is expressed or not, and the level of expression, without any negative biological significance. This constraint ensures that the analysis results have clear biological interpretability, meaning that the expression characteristics of the query cell can be regarded as a mixture of multiple homozygous cell types in a certain proportion, or approximate the expression pattern of a single cell type, avoiding meaningless negative contribution coefficients and excluding unreasonable cell type combination prediction results. By solving M NNLS models in parallel, the system quickly generates M sets of regression coefficient vectors for each query cell. Each set of coefficient vectors is analyzed, and the cell type corresponding to the element with the largest value is selected as the prediction label from the perspective of that reference subset. Finally, an initial annotation set containing M prediction labels is generated for each query cell, and each label corresponds to an independent prediction result of a reference subset.

[0024] S3: Confidence assessment (confidence assessment based on integration consistency): For each cell to be queried, count the cell type that appears most frequently in the M initial annotations and calculate the proportion of the occurrences of this type to M; when the proportion is greater than or equal to a preset threshold, the cell is classified into a high-confidence cell subgroup, otherwise it is classified into a low-confidence cell subgroup. The preset threshold is 0.7-0.9, more preferably 0.8; when the ratio is exactly equal to the threshold, it is defaulted to being classified into a high-confidence cell subpopulation; the preset threshold can be dynamically adjusted according to the batch effect strength, platform differences or species differences among reference subsets, and the dynamic adjustment of the preset threshold should be based on ensuring the purity and representativeness of the high-confidence subpopulation as the core objective.

[0025] Step S3 aims to select highly reliable samples from the initial annotation results as anchor points for subsequent adaptive calibration. The core logic involves using consistency analysis of multi-reference subset prediction results to distinguish between cells less affected by batch effects and those with significant characteristics, as well as boundary cells and disturbed cells. This provides precise targeting objects for subsequent optimization. Specifically: First, statistical analysis is performed on the M initial predicted labels obtained for each query cell. First, the label with the highest frequency among all labels (i.e., the mode label) is counted. This label represents the consensus prediction of the majority of reference subsets for this cell type. Then, the proportion of the occurrence of this mode label to the total number of reference subsets M is calculated. This proportion is defined as the consistency score (also known as the confidence score) of the cell. The higher the score, the more consistent the classification judgment of different reference subsets for this cell, and the stronger the reliability of the annotation results.

[0026] Subsequently, the system compares the consistency score with a preset confidence threshold to classify cell subpopulations. When the consistency score of a query cell is greater than or equal to the threshold, it indicates that most reference subsets have reached a consensus on the classification of the cell. The cell's expression characteristics are stable, and it is less affected by specific batch effects and technical biases, making it a reliable calibration sample. Therefore, it is identified and marked as a high-confidence cell and classified into a high-confidence cell subpopulation. Conversely, if the consistency score is lower than the threshold, it indicates that different reference subsets have significant disagreements on the classification of the cell. The cell may be in a transitional state between two or more cell types, have atypical expression characteristics, or be affected by strong batch effects and individual heterogeneity, resulting in low annotation reliability. Therefore, it is marked as a low-confidence cell and classified into a low-confidence cell subpopulation.

[0027] The preset threshold is preferably 0.8. This value has been verified by a large number of experiments and can achieve the optimal balance between the purity and recall of high-confidence cells. It ensures the annotation accuracy of high-confidence subgroups while avoiding insufficient anchor samples due to excessively high thresholds. At the same time, this threshold can be dynamically adjusted according to the signal-to-noise ratio of actual data and the degree of difference between reference subsets, thereby improving the adaptability of the method.

[0028] S4: Construction of the local optimized reference spectrum: Extract gene expression data and corresponding cell type labels from high-confidence cell subpopulations, recalculate cell type-specific pseudo-bulk expression profiles, and obtain the local optimized reference spectrum matrix; the local optimized reference spectrum matrix relies only on the high-confidence information of the query data itself, has stronger biological consistency and batch adaptability, and outliers are truncated simultaneously during the standardization process to ensure the stability and reliability of the matrix.

[0029] Step S4 breaks through the limitations of traditional annotation, which relies on external references and is difficult to adapt to local data biases, by constructing a reference spectrum that integrates external authoritative annotations and local data features. This achieves a breakthrough improvement in annotation accuracy. The system extracts the original gene expression data of all cells belonging to the high-confidence cell subpopulation from the query dataset, as well as the mode cell type label determined in the third step. This label is inherited from the authoritative annotations of the external reference dataset, ensuring the biological correctness of the cell type definition. Meanwhile, the expression data of high-confidence cells carries the unique technical characteristics and biological background of the local dataset. Thus, step S4 is the core step for achieving robust annotation in this invention.

[0030] Based on the above data pairs, the system performs pseudo-bulk expression profile calculation again, and the calculation process is consistent with the first step: the expression levels of high-confidence cells are aggregated according to cell type labels (the optimal aggregation method can be selected according to the local data characteristics), and a new matrix with the dimension of feature gene number × cell type number is constructed. Then, the same column direction standardization processing as in the first step is performed on the matrix to eliminate the slight expression bias that may exist in high-confidence cells and ensure the stability of the reference profile. The resulting matrix is ​​defined as the local optimized reference spectrum matrix, which possesses unique dual attributes: on the one hand, its column labels (i.e., cell type definitions) are fully inherited from external authoritative reference datasets, strictly adhering to validated biological classification standards to avoid cell type labeling errors; on the other hand, its numerical features (i.e., gene expression patterns) are entirely derived from the query data itself, naturally carrying the batch effect characteristics, amplification bias, sequencing technology noise patterns, and biological heterogeneity of local samples unique to the query data, perfectly adapting to the feature space of local data. This local optimized reference spectrum essentially acts as an adaptive bridge connecting the idealized external reference and the real-world query data, effectively mitigating the systematic bias between the external reference and local data.

[0031] S5: Secondary calibration annotation based on local optimized reference spectrum: Using the local optimized reference spectrum matrix obtained in step S4 as the sole benchmark, non-negative least squares analysis is performed again on all cells to be queried to obtain the cell type annotation results of the second annotation; during the secondary annotation process, a verification label is added to cells with low confidence, and subpopulation clustering, differential expression analysis or comparison with external databases can be performed to identify potential new subpopulations. After the secondary annotation in step S5 is completed, it further includes: Individual source tracing analysis is performed on low-confidence cell subpopulations, and their expression profiles can be displayed as similarity heatmaps with the local optimized reference spectrum matrix. High-confidence cell subpopulations can be further subdivided into functional subpopulations by using graph clustering or density clustering methods. The source tracing analysis further includes comparing low-confidence cells with external databases (such as CellMarker and PanglaoDB) to identify potential novel cell types or transitional cells.

[0032] Step S5 uses a localized benchmark to refine the full query data, correcting biases in the initial annotation, improving the annotation accuracy of low-confidence cells, and further verifying the annotation results of high-confidence cells, thus achieving accurate annotation of the entire dataset. The system discards the M external reference matrices used in the initial stage and uses only the local optimized reference spectrum matrix generated in step S4 as the sole benchmark. This completely eliminates the interference of bias between external references and local data, allowing the annotation process to be based entirely on a reference space adapted to local features.

[0033] Subsequently, for the complete query dataset, including high-confidence cells and low-confidence cells, an NNLS model was constructed and solved for each cell. The mathematical expression and constraints of the model were consistent with those in step S2 to ensure the coherence of the analysis logic and biological interpretability. In this process, since the baseline matrix (locally optimized reference spectrum matrix) and the vector to be solved (query cell expression vector) both came from the same sequencing batch, platform, and sample background, the systematic bias between the two was eliminated to the maximum extent, and the comparability of gene expression features was significantly improved. Therefore, the NNLS model can more accurately resolve the contribution coefficient of cell type.

[0034] This optimization is particularly effective in correcting annotations for cells with low confidence. For cells that are misclassified in the initial stage due to batch interference or atypical expression features, the local optimized reference spectrum can provide a comparison benchmark that is more in line with their actual characteristics, effectively correcting prediction bias. For transitional cells, the contribution coefficient can also accurately reflect their mixed features. Finally, the system selects the cell type with the highest contribution as the final type of each cell based on the coefficient vector obtained from the secondary analysis, while recording the confidence score of the secondary annotation and outputting the final annotation result after full calibration.

[0035] S6: Output Results: The first and second rounds of annotation results are integrated into the final annotation results, confidence scores, and local optimized reference spectrum matrix, and written to a structured file in CSV or rds format to support efficient storage and subsequent analysis; optionally, the annotation results and confidence distribution can be visualized in a low-dimensional embedding space.

[0036] Step S6 aims to provide comprehensive, traceable, and easily interpretable annotation results, taking into account the needs of subsequent data analysis and result verification, and improving the practicality of the method.

[0037] The final determined cell type labels, confidence scores of secondary annotations, M predicted labels and consistency scores of initial annotations, local optimized reference spectrum matrix generated during the construction process, and full-process operation logs (including parameter settings, preprocessing details, model solution results, and threshold selection criteria) are encapsulated and written into structured storage files, supporting multiple common formats such as CSV, rds, and h5ad. Among them, rds and h5ad formats can retain the matrix structure and meta-information of the data, which is convenient for subsequent in-depth analysis using R and Python, while CSV format is suitable for quick viewing and basic statistics.

[0038] Furthermore, to intuitively present the rationality of the annotation results and cell distribution characteristics, based on the final annotation results, dimensionality reduction algorithms such as UMAP (Uniform Manifold Approximation and Projection) or t-SNE (t-Distributed Random Neighborhood Embedding) are used to project high-dimensional gene expression data into two-dimensional or three-dimensional space, generating visual coordinates. In the graphical user interface (GUI), the system displays cell distribution in the form of a scatter plot, with different cell types marked with different colors. Confidence score information (such as represented by the size or transparency of the dots) can also be overlaid, allowing users to interactively view cell distribution, classification boundaries, and confidence distribution. In addition, auxiliary charts such as cell type composition ratio maps and confidence distribution histograms can be generated to help researchers quickly assess the rationality of the annotation results and discover potential anomalous samples or cell subpopulations.

[0039] To further improve the performance, adaptability, and computational efficiency of the above methods and meet the single-cell annotation requirements in different scenarios, this invention also includes the following steps: First, before performing step S1 to construct the initial reference spectrum, the method also includes a feature gene optimization step. The core purpose is to filter out redundant information and noise, improve model solution efficiency and annotation accuracy, and call the COSG algorithm to perform cell type-specific differential expression analysis based on all merged reference datasets. This algorithm can accurately identify the specific marker genes of each cell type, effectively avoiding the marker gene bias caused by batch effect in traditional differential analysis. It optimizes and retains genes with a fold change greater than 1.5 between different cell types, significant differences in expression rate (P value less than 0.05, corrected by multiple tests), and high expression in the target cell type, and constructs a feature gene list. This list only contains core genes with cell type differentiation ability, which greatly reduces the gene dimension of subsequent analysis.

[0040] Subsequently, based on the list of characteristic genes, a synchronous subsetting operation is performed on all reference spectrum matrices and the dataset to be queried, retaining only the genes in the list for subsequent NNLS operations. This operation not only effectively filters out the interference of housekeeping genes, background noise genes and low-difference genes, but also significantly reduces the matrix dimension, thereby improving the separability between cell types and greatly shortening the model solution time, thus improving the overall efficiency of the process.

[0041] Secondly, the nonnegative least squares (NNLS) solver has undergone deep optimization at the underlying implementation level, focusing on overcoming the computational efficiency bottleneck of large-scale single-cell datasets. The core solution algorithm is written in R language, relying on R's rich numerical computing library to ensure solution accuracy. At the same time, it combines OpenMP multithreading technology to achieve data-parallel matrix operations—distributing query cells to different CPU cores in batches to execute model construction and solving in parallel, making full use of modern multi-core CPU resources. This optimization enables the method described in this invention to have extremely strong large-scale data processing capabilities, supporting the completion of the entire annotation process within a minute-level time window for single-cell scales of more than 10,000 or even 100,000 cells. Compared with traditional serial NNLS solving, the efficiency is improved by 5-10 times, which is fully adapted to the needs of large-scale datasets generated by single-cell sequencing technology.

[0042] In addition, after the solution is completed, the obtained weight coefficient vector β is normalized (using L1 norm normalization) so that the sum of all coefficients is 1. This process gives the coefficients a clear probabilistic meaning (i.e. the contribution ratio of each cell type), which makes it easier to interpret the results. At the same time, it avoids all zero coefficients or extreme coefficients caused by extremely low gene expression levels, thus eliminating the risk of computational anomalies.

[0043] Furthermore, the confidence assessment threshold has dynamic adjustment capabilities, which can adaptively optimize according to data characteristics, improving the adaptability and robustness of the method. By calculating similarity indicators between reference subsets (such as the mean cosine similarity based on the reference spectrum matrix and the Pearson correlation coefficient), the degree of difference between reference subsets is automatically judged, and the confidence threshold is adjusted accordingly: when a small batch effect and high similarity (mean cosine similarity greater than 0.85) are detected between reference subsets, it indicates good consistency of the reference data, and the threshold can be lowered to 0.7 to improve the recall rate of high-confidence cells, obtain more anchor samples, and enhance the stability of the locally optimized reference spectrum; when significant cross-platform differences and obvious batch effects are detected (mean cosine similarity less than 0.7), it indicates poor consistency of the reference data, and the threshold needs to be raised to 0.9 to prioritize the purity of high-confidence subgroups, prevent erroneous annotation information from being introduced into the locally optimized reference spectrum, and avoid deviations in subsequent secondary annotations; when the similarity between reference subsets is at a moderate level, the default threshold of 0.8 is maintained to balance purity and recall.

[0044] Finally, for cells with extremely low confidence levels after secondary annotation (e.g., confidence scores less than 0.5), the method also includes verification labeling and source tracing analysis steps to provide auxiliary support for discovering novel cell subpopulations. These cells are uniformly labeled as Unknown or Unassigned to avoid annotation errors caused by forced classification. At the same time, source tracing analysis is automatically performed on these unknown cells: on the one hand, their expression profiles are mapped to a similarity heatmap with the locally optimized reference spectrum matrix to intuitively display the expression similarity between unknown cells and each known cell type, helping to determine their possible cell type affiliation or mixed characteristics; on the other hand, independent unsupervised clustering analysis (using the Louvain algorithm or Leiden algorithm) is performed on the subset of unknown cells, combined with differential gene analysis, to mine their unique expression characteristics, assisting researchers in discovering potential novel cell subpopulations, cells in transitional cell differentiation, doublets (i.e., two cells mistakenly captured as one), or abnormal cells affected by extreme noise.

[0045] This invention also provides a robust single-cell annotation system based on a multi-reference dataset, used to implement the aforementioned robust single-cell annotation method based on a multi-reference dataset. The system logically includes the following functional modules: Initial annotation module: configured to execute steps S1 to S3 of the above method; this module is responsible for receiving raw data, performing data preprocessing and feature gene screening; managing multiple reference subsets and scheduling the NNLS solver in parallel to perform initial global annotation; and performing statistical analysis based on voting consensus to divide the query cells into two subgroups: high confidence and low confidence.

[0046] Local reference construction module: configured to perform step S4 of the above method; this module includes a data extraction unit and a pseudo-bulk calculation unit, responsible for extracting features from the high-confidence part of the query data and dynamically generating a local optimized reference spectrum matrix in combination with cell type labels.

[0047] Secondary annotation module: configured to execute step S5 of the above method; this module is configured to load the local optimized reference spectrum matrix as a benchmark, call the NNLS solver again to perform high-precision re-annotation on the full query data, and output the final classification results and coefficient matrix.

[0048] Results visualization and interaction module: configured to execute step S6 of the above method; this module contains a dimensionality reduction algorithm library (including but not limited to UMAP, t-SNE) to map high-dimensional annotation results to low-dimensional space; and contains an interactive drawing engine that allows users to select specific regions on the visualization map and export the corresponding cell ID list or expression spectrum matrix; The results visualization module further includes a confidence level slider, customizable cell type coloring schemes, and dynamic legend generation to improve user experience and result interpretation efficiency.

[0049] Log and Audit Module: Configured to monitor and record the system's running status throughout the entire process; the content recorded by this module includes a unique task ID, the version number of the reference database used, key parameter configurations (such as confidence threshold, number of feature genes), the total process time, and the hardware and software environment information during runtime.

[0050] It also includes a user permission management module to support concurrent access by multiple users and data isolation.

[0051] A specific application embodiment of the present invention is as follows: This embodiment demonstrates the application of the robust single-cell annotation method (RoMCA) described in this invention to a publicly available human peripheral blood mononuclear cell (PBMC) dataset and compares its performance with existing technologies.

[0052] 1. Experimental Environment and Data Preparation The algorithm implementation in this embodiment is based on the R language environment.

[0053] First, load the necessary R packages, including Seurat, nnls, COSG, and other auxiliary packages.

[0054] The data preparation process is as follows: (1) Read in the wilk Seurat object, which contains PBMC single-cell data from multiple donors, and each cell has a real cell type label; (2) Perform normalization (NormalizeData) and high-variance gene (FindVariableFeatures) calculations on the wild objects; (3) Split the wild object according to Donor (sample source) to obtain wild_list; (4) Specify a donor as the query dataset to be annotated; (5) Subset the data from the remaining donors or keep it in list form as a reference dataset (wilk_ref or wild_list). (6) Extract the true labels of the wild_query for subsequent accuracy evaluation.

[0055] 2. Implementing the RoMCA annotation method of the present invention The core process of this invention (such as) Figure 1 (As shown) Follow these steps: Step 1: Initial Reference Spectrum Construction and Initial Annotation (1) Constructing a multi-reference spectrum: Iterate through the reference dataset list `wilk_list` (i.e., the reference set split by Donor). In the loop, for each Donor reference subset: (2) Obtain its expression matrix ref_expr and metadata ref_meta: By calculating the mean gene expression of each cell type, construct the pseudo-bulk reference spectrum ref_mean_mat of the Donor. Gene screening functions (such as get_marker_gene and second_filter) can be called to optimize the reference spectrum to enhance its specificity.

[0056] (3) Initial annotation: Obtain the expression matrix query_expr of the query dataset wilk_query. In the above loop, call the non-negative least squares annotation function nnls_prd defined in this invention. The nnls_prd function accepts the reference spectrum ref_mean_use (as matrix A) and the query data query_use (as matrix B). The function solves by nnls(A, B) and returns the cell type corresponding to the maximum coefficient value after normalization. Then, the annotation results of each Donor reference spectrum for the query set are stored in the list prd_list.

[0057] Step 2: Confidence Assessment (1) After the initial annotation loop ends, prd_list is converted into a data frame prd_df. At this time, each row of prd_df represents a query cell and each column represents an annotation result from a different Donor reference spectrum.

[0058] (2) Call the get_mode function to count the most frequent annotation results in each row (each cell), and use them as the initial high-confidence prediction prd_df. anno1.

[0059] (3) Calculate the highest frequency value using get_mode2 (or a similar method).

[0060] (4) Set a confidence threshold (e.g., 0.8) and mark the annotations (anno2) of cells with a frequency below the threshold as unknown.

[0061] Step 3: Local Reference Building (1) such as Figure 1 As shown in the flowchart, this invention filters out all cells that were not marked as unknown in step two from the query dataset wilk_query, i.e., the high-confidence cell subpopulation.

[0062] (2) Extract the expression data ref_expr of these high-confidence cells and their reliable annotation ref_meta.

[0063] (3) Using these cells derived from the query data itself, repeat the pseudo-batch construction method in step three to generate a new, highly adapted local optimized reference spectrum matrix.

[0064] Step 4: Secondary Annotation (1) such as Figure 1 As shown in the flowchart, obtain the expression matrix query_expr of the complete query dataset wilk_query.

[0065] (2) Call the nnls_prd function again. This time, the reference spectrum (matrix A) is the local reference spectrum constructed in step three, and the query data is the complete query set.

[0066] (3) The obtained pred_celltype is the secondary annotation result. It is stored in prd_df[,'anno2'] and integrated with the first round of annotation results to output the final annotation result of the present invention.

[0067] 3. Performance Evaluation and Comparison (1) Results of the RoMCA method of the present invention Compare the final annotation results obtained in step four with the actual labels.

[0068] On the H2 sample (i.e., wild_query) in this embodiment, the annotation accuracy of the RoMCA method of the present invention reached 93.2%.

[0069] like Figure 2 As shown in the figure, the confusion matrix heatmap of the RoMCA method of the present invention shows that the prediction results are highly concentrated on the diagonal, indicating that all cell types are accurately identified and that there is extremely high discrimination between easily confused cell subtypes such as CD14 Monocyte, CD16 Monocyte, CD4n T, and CD8m T.

[0070] (2) Comparison with other algorithms To verify the beneficial effects of the present invention, this embodiment reproduces a mainstream cell type algorithm, SciBet, on the exact same query set and reference set; on the same H2 samples, the annotation accuracy of the existing SciBet algorithm is 90.7%.

[0071] Comparative results show that the accuracy of existing technologies is limited by the technical differences between the reference and query sets and by batch effects; the RoMCA method of this invention (accuracy of 93.2%) is significantly better than existing technologies, such as SciBet's 90.7%; therefore, the beneficial effect of this invention is that it innovatively adopts a two-step strategy of local reference construction and secondary annotation, and uses the high confidence information of the query data itself to dynamically calibrate the annotation model, effectively overcoming the batch effect, thereby achieving more robust and accurate cell type annotation.

[0072] To comprehensively and preliminarily verify the robustness and importance of the RoMCA algorithm proposed in this invention in handling batch effects in single-cell RNA sequencing (scRNA-seq) data, this embodiment designed rigorous cross-dataset benchmark tests. In this test, Seuratv4 (a point integration-based label transfer algorithm) and SciBet (an entropy-based fast classification algorithm), both recognized mainstream algorithms in the field of bioinformatics, were tested multiple times as benchmark products.

[0073] Dataset and Experiment Setup Data source: This embodiment uses two human peripheral blood mononuclear cell (PBMC) datasets containing multiple donor sources, one as a reference dataset and the other as a query dataset; both datasets have significant inter-sample heterogeneity and are ideal benchmarks for evaluating the robustness of the algorithm.

[0074] Evaluation metric: Accuracy, which is the proportion of all correctly predicted samples out of the total sample.

[0075] Comparison Algorithm Three cell type annotation algorithms, SciBet, Seurat v4, and singleR, were used. To ensure fairness in the comparison, all algorithms were run with default parameters.

[0076] Test results are available Figure 3 The specific data quantification is as follows: RoMCA: Accuracy rate is 99.24%.

[0077] Scibet: Accuracy rate is 96.47%.

[0078] SingleR: Accuracy rate was 69.16%.

[0079] Seurat: Accuracy rate 98.96% according to Figure 3 Experimental results show that the RoMCA algorithm of this invention achieves the best performance in cross-dataset testing, outperforming existing mainstream tools.

[0080] Based on the benchmark test results of the above comprehensive dataset, the RoMCA algorithm of this invention shows an accuracy advantage over mainstream cell type annotation algorithms. RoMCA actually uses high-quality information within the query set to guide the classification of low-quality cells, significantly improving the overall data utilization and annotation completeness.

[0081] In summary, cross-dataset benchmark tests fully demonstrate that the RoMCA algorithm significantly outperforms current industry standard algorithms in terms of accuracy, robustness, and understanding of complex biological samples.

[0082] This implementation of the RoMCA algorithm can be deployed on any electronic device with computing power (including but not limited to servers, workstations, or personal computers) to automate the single-cell annotation process. The algorithm includes three core functional modules, each implemented through software code (such as R scripts) and working collaboratively under the control of the computer processor: Initial annotation module: This module is mainly responsible for data preprocessing and preliminary analysis. Its specific functions include: Multi-reference management: Receives user-inputted reference datasets and splits them into multiple independent reference subsets based on sample source; Pseudo-batch construction: For each reference subset, calculate the average expression profile of each cell type and construct the initial reference matrix; Parallel NNLS parsing: Receives the single-cell expression matrix to be queried and calls the non-negative least squares (NNLS) algorithm in parallel to calculate the cell type coefficient of each queried cell in each reference subset and outputs a preliminary annotation result list.

[0083] Secondary annotation module: designed to eliminate batch effects; this module automatically extracts high-confidence cells and their corresponding high-confidence labels from the query dataset, which have been filtered by the previous module. Using these extracted data, this module dynamically calculates a new pseudo-batch expression profile, namely the local optimized reference profile, which is biologically homologous to the query data, thereby minimizing technical noise.

[0084] Results visualization and interaction module: This module summarizes the results of the first two rounds of annotation, and considers cell types with the same annotation results in the initial annotation module and the secondary annotation module as reliable cell type annotation results and outputs them.

[0085] Based on the above system architecture, this embodiment can automatically complete the entire process from data import to final annotation when processing PBMC data containing complex batch effects. Furthermore, each module is not redundant during cell type annotation. Ablation experiments were performed on the system described in this embodiment, and the results are as follows: Figure 4 As shown, through Figure 4 The results of the ablation experiments showed that the model with the three modules stacked together had a higher annotation accuracy than the single fine module.

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

Claims

1. A robust single-cell annotation method based on a multi-reference dataset, characterized in that, Includes the following steps: S1: Initial reference spectrum construction: Obtain single-cell RNA sequencing reference datasets from two or more independent biological samples, and label them with cell type tags. Divide the reference datasets into M reference subsets (M≥3). For each reference subset, group them according to the labeled cell type tags, calculate the average gene expression value of each group of cells, and obtain the pseudo-bulk expression matrix of cell type × gene, which is denoted as the initial reference spectrum matrix. S2: Initial annotation: Extract the gene expression vector of each cell in the single-cell sequencing dataset to be queried. Based on the M initial reference spectrum matrices obtained in step S1, solve the non-negative least squares model to obtain the M initial cell type annotations for each cell. S3: Confidence assessment: For each cell to be queried, count the cell type that appears most frequently in the M initial annotations and calculate the proportion of the number of times that type appears in M. When the proportion is greater than or equal to a preset threshold, the cell is classified into a high-confidence cell subpopulation; otherwise, it is classified into a low-confidence cell subpopulation. S4: Construction of the local optimized reference spectrum: Extract gene expression data and corresponding cell type labels from high-confidence cell subpopulations, recalculate cell type-specific pseudo-bulk expression profiles, and obtain the local optimized reference spectrum matrix; The locally optimized reference spectrum matrix relies solely on the high-confidence information of the query data itself, exhibiting stronger biological consistency and batch adaptability. Outliers are truncated during the standardization process to ensure the stability and reliability of the matrix. S5: Secondary calibration annotation based on local optimized reference spectrum: Using the local optimized reference spectrum matrix obtained in step S4 as the sole benchmark, non-negative least squares analysis is performed again on all cells to be queried to obtain the cell type annotation results of the second annotation. During the secondary annotation process, cells with low confidence are attached with a marker to be verified, and subpopulation clustering, differential expression analysis, or comparison with external databases can be performed to identify potential new subpopulations. S6: Output Results: The first and second rounds of annotation results are integrated into the final annotation results, confidence scores, and local optimized reference spectrum matrix, and written to a structured file in CSV or rds format to support efficient storage and subsequent analysis; optionally, the annotation results and confidence distribution can be visualized in a low-dimensional embedding space.

2. The robust single-cell annotation method based on a multi-reference dataset according to claim 1, characterized in that, The splitting method must ensure the independence and diversity of the reference subsets. The reference subsets are split according to biological individuals, sequencing batches, or platform sources, and the number of cells in each subset is ≥500. Furthermore, for subsets with fewer than 500 cells, random oversampling or merging into adjacent subsets is used to ensure statistical robustness. The splitting method further includes dividing according to disease state, tissue location, or time point to increase the biological diversity of the reference subsets.

3. The robust single-cell annotation method based on a multi-reference dataset according to claim 1, characterized in that, The objective function of the nonnegative least squares model is min‖y-Xβ‖. 2 2. The constraint β≥0, where y is the gene expression vector of a single query cell, X is the reference spectrum matrix, and β is the cell type contribution coefficient vector. After solving the problem, the cell type corresponding to max(β) is used as the annotation result under this reference subset. Furthermore, β is normalized; The model is solved using a low-level R implementation to support efficient parallel computation, accelerated by pbmclapply, and supports computations > 10. 4 Rapid processing of large numbers of cells.

4. The robust single-cell annotation method based on a multi-reference dataset according to claim 1, characterized in that, The preset threshold is 0.7-0.9, more preferably 0.8; when the ratio is exactly equal to the threshold, it is defaulted to being classified into a high-confidence cell subpopulation; the preset threshold can be dynamically adjusted according to the batch effect strength, platform differences or species differences among reference subsets, and the dynamic adjustment of the preset threshold should be based on ensuring the purity and representativeness of the high-confidence subpopulation as the core objective.

5. A robust single-cell annotation method based on a multi-reference dataset according to claim 1, characterized in that, After the secondary annotation in step S5 is completed, it further includes: Individual source tracing analysis is performed on low-confidence cell subpopulations, and their expression profiles can be displayed as similarity heatmaps with the local optimized reference spectrum matrix. High-confidence cell subpopulations can be further subdivided into functional subpopulations by using graph clustering or density clustering methods. The source tracing analysis further includes comparing low-confidence cells with external databases to identify potential novel cell types or transitional cells.

6. The robust single-cell annotation method based on a multi-reference dataset according to claim 1, characterized in that, The gene screening step is included before step S1: (1) The COSG differential expression algorithm was used to identify cell type-specific marker genes from the merged reference dataset; (2) Perform fold change analysis on the pseudobulk expression profile and retain genes with FoldChange>1.5 among different cell types; (3) Using the gene list obtained in steps S1 and S2, the reference spectrum matrix and the dataset to be queried are synchronized and subsetted to reduce noise and improve computational efficiency. The gene screening step further includes removing genes with low expression or high deletion rates to ensure the quality and robustness of the input data.

7. A robust single-cell annotation system based on a multi-reference dataset, used to implement the robust single-cell annotation method based on a multi-reference dataset as described in any one of claims 1-6, characterized in that, include: The initial annotation module is used to perform steps S1 to S3 in claim 1 and output high / low confidence cell subpopulations. The local reference building module is used to execute step S4 and generate a local optimized reference spectrum matrix; The secondary annotation module is used to execute step S5 and output the final annotation result; The results visualization module is used to display annotation results and confidence distribution in a low-dimensional space; The logging and auditing module records the parameters, version number, and execution time of each annotation task to ensure that the results are traceable. It also includes a user permission management module to support concurrent access by multiple users and data isolation.

8. A robust single-cell annotation system based on a multi-reference dataset according to claim 7, characterized in that, Both the initial annotation module and the secondary annotation module are implemented using the same non-negative least squares solver, and the solver supports parallel computation to accelerate processing >10. 4 Processing of cells at the scale.

9. A robust single-cell annotation system based on a multi-reference dataset according to claim 7, characterized in that, The results visualization module supports multiple dimensionality reduction algorithms, including UMAP and t-SNE, and supports interactive selection and export of high-confidence cell subpopulations. The results visualization module further includes a confidence level slider, customizable cell type coloring scheme, and dynamic legend generation to improve user experience and result interpretation efficiency.

Images