Single cell data quality control processing method and system

Abstract

The invention relates to a single cell data quality control processing method. The method comprises the following steps: 1) calculating a quality control index according to single cell data to be processed; identifying outliers in the single cell data to be processed according to the quality control indexes and / or the cell type composition to obtain a sample level quality control analysis result; 2) identifying empty liquid drops and double-cell data in the single-cell data to be processed; then identifying low-quality cell data in the to-be-processed single cell data to obtain a cell level quality control analysis result; 3) calculating and evaluating the average expression level, the expression standardization variance and the detection rate in all cells of each feature in the to-be-processed single cell data to obtain a feature level quality control analysis result; and 4) evaluating the batch effect of the to-be-processed single cell data to obtain a batch effect evaluation analysis result. Compared with other existing process tools, the single-cell multi-omics quality control process tool is more complete, more flexible and more friendly, and provides the most comprehensive and powerful quality control workflow for single-cell multi-omics analysis.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a method and system for quality control processing of single cell data. Background Art

[0002] The development of single-cell omics technology has played a key role in promoting people's understanding of health and disease, and has provided an important basis for discovering disease mechanisms and potential therapeutic targets. Quality control (QC) of single-cell multi-omics is crucial to ensure the reliability and accuracy of data. In order to ensure data integrity, avoid mislabeling and misidentification, and thus prevent misleading results in downstream analysis, it is crucial to systematically evaluate low-quality items in samples, cells, and expression features, as well as batch effects from various technologies and human errors, environmental contamination, and background noise.

[0003] Although there have been significant advances in single-cell QC strategies, for example, scDblFinder and DoubletFinder can reliably identify and exclude double cells; tools such as CellBender and DecontX can effectively alleviate system background noise; SCTK-QC workflow can generate and visualize QC indicators for scRNA-seq data. However, existing quality control tools and methods have the following limitations: First, existing technologies are mainly targeted at a single group of omics (such as single-cell transcriptomics) and lack comprehensive quality control of multiple omics. In addition, the quality control process is relatively fragmented and lacks a systematic and comprehensive quality control evaluation. In addition, there is a lack of quality control methods for comparison between samples, which are usually limited to single sample quality control.

[0004] Therefore, in view of the defects of current quality control methods, the present invention proposes a complete and standardized single-cell data quality control processing method. Summary of the invention

[0005] In order to perform systematic and comprehensive quality control of single-cell data in multi-omics and ensure the consistency and biological accuracy of multi-dimensional data, the present invention proposes a method for quality control processing of single-cell data. The present invention is implemented by the following technical solutions:

[0006] A method for single cell data quality control processing comprises the following steps:

[0007] Step 1) calculating quality control indicators in the single-cell data to be processed; identifying outliers in the single-cell data to be processed based on the quality control indicators and / or cell type composition, and obtaining sample-level quality control analysis results;

[0008] Step 2) identifying empty droplets and double cell data in the single cell data to be processed; then identifying low-quality cell data in the single cell data to be processed to obtain cell-level quality control analysis results;

[0009] Step 3) Calculate the average expression level, expression normalization variance and detection rate of each gene / protein in all cells in the single-cell data to be processed, and evaluate and remove the background noise of the single-cell data to obtain feature-level quality control analysis results;

[0010] Step 4) Evaluate the batch effect of the single-cell data to be processed and obtain batch effect evaluation analysis results.

[0011] The method of the present invention can process single-cell single-omics data as well as single-cell multi-omics data.

[0012] Optionally, the single-cell multi-omics includes: scRNA-seq, RNA+ADT-seq and scTCR / BCR-seq.

[0013] Preferably, the quality control indicators in step 1) include: cell number, median UMI abundance of sample cells, median gene / protein detected in sample cells, median mitochondrial gene count percentage of sample cells, TCR% and BCR%, etc.

[0014] Optionally, the step 1) identifying outliers in the single cell data to be processed based on the quality control index comprises the following steps:

[0015] The median absolute deviation is used to calculate the deviation of the value of each quality control indicator in the processed single-cell data from the median, and the outlier samples based on each quality control indicator are identified;

[0016] Preferably, the quality control indicators include the median of UMI abundance in sample cells, the median of detected genes / proteins in sample cells, and the median of mitochondrial gene count percentage in sample cells.

[0017] Optionally, in step 1), identifying outliers in the single-cell data to be processed based on cell type composition comprises the following steps:

[0018] 1) Compare and analyze the cell type ratio of the sample in the single-cell data to be processed with the cell type ratio of the internally constructed tissue atlas to determine whether the cell ratio of the sample deviates from the overall range of its cell type, and identify the deviated samples as outliers;

[0019] 2) PCA combined with confidence intervals and DBSCAN method was used to identify samples with outliers in the overall cell type ratio of the sample group.

[0020] Optionally, the step 1) further includes performing identity authentication processing on the single-cell data to be processed, and identifying samples with identity authentication errors;

[0021] Preferably, the identity verification process includes setting a verification gene based on the sample type of the single-cell data to be processed, and performing identity verification by evaluating the cell detection ratio of the verification gene.

[0022] Optionally, in step 2), the dropletUtils package is used to identify empty droplets;

[0023] Preferably, in the step 2), at least one of a single-cell transcriptome dual-cell identification method, a V(D)J sequencing dual-cell identification method, and an ADT data dual-cell identification method is used to perform dual-cell data identification;

[0024] Preferably, the single-cell transcriptome dual-cell identification method includes at least one of DoubletFinder, scDbFinder, cxds, bcds, and hybrid;

[0025] Preferably, the ADT data dual-cell identification method is to perform dual-cell identification based on the ADT data in the RNA+ADT combined sequencing technology by combining the cyclic iterative COSG method and the common statistical classification method.

[0026] Optionally, the step of identifying low-quality cell data in the single-cell data to be processed adopts at least one of fixed threshold filtering, MAD filtering, miQC filtering, and ddqc filtering.

[0027] Optionally, in the step 2), the identified low-quality data are further subjected to cluster analysis, and the cell types of the identified low-quality cell data are analyzed according to the cluster analysis results;

[0028] Then determine whether to retain this part of the single cell data based on the cell type.

[0029] The feature in step 3) is the gene or protein data in the single-cell data, and each gene and protein is a feature.

[0030] Optionally, the step 3) further includes performing background noise correction on the single cell data to be processed, and the background noise correction includes: global background noise correction or local background noise correction;

[0031] Preferably, the global background noise correction method in step 3) includes DecontX and DecontPro methods to achieve global background noise correction for RNA and ADT data.

[0032] Preferably, the local background noise correction method in step 3) includes using scCDC to identify and remove highly contaminated RNA or protein.

[0033] Optionally, the evaluating the batch effect of the single-cell data to be processed includes:

[0034] 1) Perform UMAP analysis on the processed single-cell data and PCA analysis on the pseudo-transcriptome data;

[0035] 2) Singular value decomposition analysis and scater package are used to analyze the influence of batch-related variables.

[0036] Optionally, the method for single cell data quality control processing further includes:

[0037] The sample-level quality control analysis results, cell-level quality control analysis results, feature-level quality control analysis results and batch effect evaluation analysis results are combined into a visual web report.

[0038] A method for single cell data quality control processing, comprising:

[0039] A sample-level quality control module is used to calculate quality control indicators in the single-cell data to be processed; identify outliers in the single-cell data to be processed based on the quality control indicators and / or cell type composition, and obtain sample-level quality control analysis results;

[0040] The cell-level quality control module is used to identify empty droplets and double cell data in the single-cell data to be processed; then identify low-quality cell data in the single-cell data to be processed, draw quality control charts, and obtain cell-level quality control analysis results;

[0041] The feature-level quality control module is used to calculate and evaluate the average expression level, expression normalization variance and detection rate of each gene / protein in the single-cell data to be processed, and to evaluate and remove the background noise of the single-cell data to obtain the feature-level quality control analysis results;

[0042] The batch effect analysis module is used to evaluate the batch effect of the single-cell data to be processed and obtain the batch effect evaluation analysis results.

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

[0044] The present invention can simultaneously perform quality control on multi-omics sequencing data of multiple samples, including scRNA-seq, surface protein sequencing, and immune repertoire sequencing, allowing users to fully evaluate their data before downstream analysis. The present invention integrates existing quality control steps, builds a bridge between various quality control processes and tools, solves the compatibility issues between different tools, and connects multiple quality control processes into one workflow. The present invention is built on the most popular Seurat framework, which is convenient for users to perform downstream analysis. At the sample level, the present invention provides a unique outlier sample detection module, especially a cell ratio outlier detection module, which assists users in judging the biological variation of samples by integrating thousands of single-cell data, and assists in sample-level quality control. At the cell level, the present invention provides low-quality cell filtering standards and recommended reference ranges for sequencing features for various tissues, enabling users to select the most appropriate QC reference threshold for their experiments. In addition, in order to avoid excessive filtering of non-low-quality cells, the present invention also provides a re-clustering strategy for low-quality cells, as well as marker gene enrichment analysis. In addition, the present invention also develops a dual-cell identification strategy for ADT data in the multi-omics RNA+ADT joint sequencing technology to make up for the shortcomings of dual-cell identification at the RNA level.

[0045] The present invention provides a large number of integrated and convenient calculations and QC charts, which greatly simplifies the complex quality control process. In particular, the present invention can generate interactive HTML reports containing single-cell data system quality control information, allowing users to quickly check data quality before downstream analysis, allowing users who do not need in-depth programming experiments to easily perform comprehensive data quality control.

[0046] The present invention is compatible with Seurat V5 and BPCells versions, can adapt to large-scale single-cell data quality control, and reduce memory usage. It has low requirements for computer memory, and an ordinary laptop can perform quality control analysis on millions of cells.

[0047] Compared with other existing process tools, the single-cell multi-omics quality control process tool of the present invention is more flexible and user-friendly, providing the most comprehensive and powerful quality control workflow for single-cell multi-omics analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0049] Figure 1 It is a flow chart of the present invention.

[0050] Figure 2 It is the basic situation of the data, where A is the quality control index of the data; B is the proportional distribution of TCR and BCR chains in V(D)J data; and C is the receptor situation distribution of V(D)J data.

[0051] Figure 3 It is the outlier sample identification situation, where A is the outlier situation of the basic quality control indicators of the sample; B is the violin distribution diagram of the basic quality control indicators of the sample; C is the PCA distribution diagram of the proportion of cell types in the sample.

[0052] Figure 4 is the result of two-cell recognition, where A is the proportion of two-cells identified by different two-cell recognition methods; B is the similarities and differences of the prediction results of different two-cell methods.

[0053] Figure 5 It is the proportion of low-quality cells identified by different low-quality cell identification methods.

[0054] Figure 6 It is the result of low-quality cell re-clustering analysis, where A is the result of automatic annotation of low-quality cell re-clustering ScType; B is the situation of some classic cell markers in different re-clustering cell groups.

[0055] Figure 7 It is the result of feature hierarchical analysis, where A is the expression of cell genes in the cell population; B is the expression of highly polluted genes identified by the scCDC method in the cell population; C is the expression of cell proteins in the cell population; and D is the expression of highly polluted proteins identified by the scCDC method in the cell population.

[0056] Figure 8 These are partial analysis results of batch effects, where A is the PCA analysis at the sample level; B is the UMAP map of RNAomics at the cell level; and C is the UMAP map of ADTomics at the cell level.

[0057] Fig. 9 is the result of batch covariate analysis, where A is the SVD analysis of batch covariates and PC1 of pseudo transcriptome data; B is the contribution of each batch covariate at the PC1 level to the pseudo transcriptome data;

[0058] Fig.10 This is an example of a report page generated by the system of the present invention.

[0059] Fig.11Figure 3 is a scatter plot of protein expression threshold division without background correction, where A is the scatter plot of expression threshold division of mutually exclusive protein pairs of CD3-CD20 and CD19-CD3 of TP1 sample, and the three figures correspond to different method strategies. Red represents the division threshold of each method, of which the first two are the strategies of this method, and the last one is the MLTiplet CITE-seq method strategy. The red dot ambiguous represents the TB double cells identified at the VDJ level; B is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP2 samples; C is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP3 samples; D is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP3-rep samples; E is the scatter plot of expression threshold division of mutually exclusive protein pairs of pbmc5K samples; and F is the scatter plot of expression threshold division of mutually exclusive protein pairs of pbmc10K samples.

[0060] Fig.12 It is a scatter plot of protein expression threshold division under background correction; A is the scatter plot of expression threshold division of CD3-CD20 and CD19-CD3 mutually exclusive protein pairs of TP1 sample, and the two figures correspond to different method strategies, red represents the division threshold of each method, and the red dot ambiguous represents the TB double cells identified at the VDJ level; B is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP2 samples; C is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP3 samples; D is the scatter plot of expression threshold division of mutually exclusive protein pairs of TP3-rep samples; E is the scatter plot of expression threshold division of mutually exclusive protein pairs of pbmc5K samples; F is the scatter plot of expression threshold division of mutually exclusive protein pairs of pbmc10K samples.

[0061] Fig.13 The threshold changes of various methods with different clustering resolutions; A is the threshold fluctuation of different methods for TP1 sample division; B is the threshold fluctuation of different methods for TP2 sample division; C is the threshold fluctuation of different methods for TP3 sample division; D is the threshold fluctuation of different methods for TP3-rep sample division; E is the threshold fluctuation of different methods for pbmc5K sample division; F is the threshold fluctuation of different methods for pbmc10K sample division. DETAILED DESCRIPTION

[0062] In order to make the purpose and technical solution of the present invention more clear, the present invention is further described in detail. The experimental methods described in the following examples are all conventional methods unless otherwise specified: if no specific techniques or conditions are specified in the examples, they are carried out according to the techniques or conditions described in the literature in the field or according to the product instructions; the reagents and materials, unless otherwise specified, can be obtained from commercial channels.

[0063] The present invention proposes a method for single cell data quality control processing, comprising the following steps:

[0064] Step 1) Sample-level quality control, including: calculating quality control indicators in the single-cell data to be processed; identifying outliers based on the quality control indicators and cell types in the single-cell data to be processed. The specific process includes:

[0065] 1-1) Sample quality assessment:

[0066] This step is used to identify samples that may be abnormal due to library construction, insufficient sequencing depth, amplification bias, and contamination. By integrating sequencing indicators generated by upstream software such as CellRanger, such as Q30, saturation, etc., the quality of the samples can be quickly and preliminarily evaluated according to the official standards. Then, other common quality control indicators are calculated. Detailed sample-level quality control indicators are shown in Appendix 1. Some of the more important ones, such as cell number (nCell), median of nCount_RNA / nCount_ADT (median of sample cell UMI abundance), median of nFeature_RNA / nFeature_ADT (median of sample cell detected genes / proteins), median of MT% (median of sample cell mitochondrial gene count percentage), TCR% and BCR%, etc., can be combined with quality control visualization charts to perform preliminary quality control of samples.

[0067] Table 1 shows all the optional quality control indicators in the method of the present invention.

[0068] Table 1: Sample-level quality control indicators

[0069]

[0070]

[0071]

[0072] 1-2) Outlier sample detection:

[0073] This step identifies samples with significant outliers in quality control indicators or cell type composition within the sample group. Outlier sample detection includes two aspects: a) and b):

[0074] a) Identify samples within the sample group that deviate from the general distribution of various quality control indicators. The assumption of this strategy is that the values ​​of sample quality control indicators are generally similar. The default quality control indicators mainly include median of nCount_RNA / nCount_ADT, median of nFeature_RNA / nFeature_ADT and median of MT% common quality control indicators. Other quality control indicators can also be selected, see Table 1. The method of identifying outliers is to use the median absolute deviation (MAD) to measure the degree of deviation of each value in the data set from the median, thereby identifying outlier samples for each key quality control indicator. Each quality control indicator will warn of potential outlier samples. The more quality control indicator warnings a sample has, the more likely it is that the sample is abnormal. This step allows users to quickly view potential abnormal samples, assist users in sample control, and avoid the impact of low-quality samples on biological analysis conclusions.

[0075] b) Identify samples with significant outliers in cell type composition. Deviations in cell type composition among samples from the same experimental group may indicate underlying technical issues such as extended processing time, changes in storage conditions, experimental processing errors, or batch effects, all of which may compromise the integrity of the data. The strategy for identifying sample cell type composition consists of two aspects:

[0076] First, the cell type ratio of the sample was compared with the cell type ratio of the existing tissue atlas. The specific steps include: first downloading 2299 single-cell samples from DISCO data, involving 10 tissue types, including blood, bone, bone_marrow, thymus, intestine, kidney, liver, lung, pancreas and skin. Each sample was clustered and annotated using the ScType automatic annotation method and the self-defined main cell type population, and the fluctuation of the main types in different tissue samples was further counted to establish the reference fluctuation range of the single cell type of each tissue. The cell types include B cells (B cell), dendritic cells (DC), endothelial cells (Endothelial), fibroblasts (Fibroblast), granulocytes (Granulocyte), monocytes / macrophages (Mon / Mac), natural killer cells (NK), other cells (Other) and T cells (T cell). The reference fluctuation range is the maximum and minimum values ​​of the cell type in all samples of the tissue (Table 2). The reference fluctuation range is the maximum and minimum values ​​of the cell type in all samples of the tissue. For each sample to be tested, clustering annotation will be automatically run and the cell type ratio will be further counted. It will be compared with the established tissue fluctuation range to finally determine whether the cell ratio of the sample deviates from the overall range. If it deviates, the sample will be identified as an outlier. The sample of the outlier is more likely to have abnormal conditions, and the user will be warned.

[0077] Table 2: Reference fluctuation ranges for single cell types in different tissues

[0078]

[0079]

[0080]

[0081] Second, use PCA combined with confidence intervals and DBSCAN methods to identify samples with outliers in the overall cell type ratio of the sample group. The assumption of this type of strategy is that samples in the same group have more similar cell ratio compositions. Specifically, based on the existing annotation results based on ScType, the cell type ratio of each sample is calculated, and PCA analysis is performed on all samples, using confidence ellipses and DBSCAN to identify potential outliers. These outlier samples represent a shift in cell ratios, which to a certain extent represents the heterogeneity of such samples, and has different meanings in different groups. Users can conduct specific analysis based on the circumstances of their respective data sets. For example, for samples from different batches of the same person, theoretically the ratio should remain consistent. If significant outliers are found, it indicates that the sample is abnormal.

[0082] 1-3) Sample Authentication Module:

[0083] Purpose: To identify incorrect sample labels that may be caused by sample handling errors.

[0084] Solution: Assess the expression pattern of specific genes by evaluating the proportion of cells expressing validated genes, and compare the consistency of data expression performance with recorded labels through visual charts. By default, users can determine the sex of samples based on sex-specific gene expression, including XIST, DDX3Y, UTY, and RPS4Y1, to better check consistency with recorded sex.

[0085] Step 2) Cell-level quality control, including identifying empty droplets and double cell data in the single cell data to be processed; then identifying low-quality cell data in the single cell data to be processed. The specific method includes the following steps:

[0086] 2-1) Identification and removal of empty droplets and double cells. It avoids result errors caused by environmental RNA contamination, microfluidic cell capture errors and improper cell encapsulation. Among them, empty droplets are identified using the dropletUtils package. The identification of double cells integrates five existing public single-cell transcriptome methods (DoubletFinder, scDbFinder, cxds, bcds and hybrid), a public method using V(D)J data, and a new double cell identification method for ADT data in RNA+ADT combined sequencing technology. Generally, the double cell detection strategy in V(D)J sequencing will mark cells with multiple TCR or BCR chains (or both chains) as potential double cells. This is particularly effective in identifying homotypic double cells by detecting multi-chain cells, which can make up for the shortcomings of transcriptome methods in identifying homotypic double cells. For single-cell data, single-cell transcriptome strategies and other omics methods provided by this process can be combined to jointly identify double cells.

[0087] Single-cell transcriptome analysis is prone to information loss due to the sparsity of data, resulting in blind spots in quality control. The ADT data in the RNA+ADT combined sequencing technology is bimodal, and the missing rate is not so serious. Therefore, for the multi-omics RNA+ADT combined sequencing technology data, the identification of dual cells can be further identified and supplemented from the protein level. The dual-cell identification method developed internally based on the ADT data in the RNA+ADT combined sequencing technology of the present invention mainly includes two major parts. The first step is to identify the negative and positive cell populations using the cyclic iterative COSG strategy, and the second step is to use statistical methods to identify the threshold level of distinguishing high and low expression at the protein level. Specifically, the dual-cell identification method comprises the following steps:

[0088] a) Obtain a data set to be subjected to dual-cell identification, and then cluster the data set to obtain multiple cell populations. The data set is a single-cell data set generated by RNA+ADT combined sequencing technology. Clustering adopts RNA clustering or ADT clustering. The present invention supports RNA and ADT single-omics clustering. When the number of ADT-omics proteins is relatively small, RNA clustering can be used instead, and ADT clustering will be used by default.

[0089] b) using the disclosed COSG method with cyclic iteration to determine the negative expression cell population and the positive expression cell population corresponding to each protein in the mutually exclusive protein pair from the cell population obtained in step a); the mutually exclusive protein pair is customized according to the biological cell type, for example, the mutually exclusive protein pair of B and T cells is CD19 and CD3; specifically comprising:

[0090] For each protein in the mutually exclusive protein pair, the percentage of the protein detected in each cell population and the average expression level are calculated; the protein percentage and the average expression level are multiplied to obtain the penalty mean of the cell population;

[0091] The cell population with the highest penalty mean is used as the initial seed, and it is iterated and added in descending order according to the penalty mean of the cell population; in each iteration, the COSG score of the combination of the added cell populations is calculated using the public COSG method. This public COSG method calculates the cosine similarity between the expression pattern of a feature in a certain cell population and its assumed marker pattern, thereby obtaining the specific COSG score of the feature in a certain cell population. The specific formula is as follows:

[0092]

[0093] Among them, g i Represents each protein, G k represents the combination of cell populations, λ k Represents G kExpression pattern of putative marker genes for cell population combinations, λ t Represents except G k The expression pattern of the assumed marker gene of other cell populations outside the COSG, u represents the penalty factor; when the COSG score does not increase after iteration or the frequency of the corresponding marker gene in the current cell population is less than 0.2, the iteration is stopped; the cell population in the added cell population combination is the expression-positive cell population, and the remaining cell populations are the expression-negative cell populations.

[0094] c) calculating the demarcation threshold between positive expression and negative expression of each protein in the mutually exclusive protein pair according to the data in the expression-negative cell population and the expression-positive cell population;

[0095] d) identifying the dual-cell data in the data set to be subjected to dual-cell identification according to the demarcation threshold. Generally speaking, the cell data whose expression level of each protein in the mutually exclusive protein pair is higher than the respective threshold is identified as dual-cell data.

[0096] d) Calculating the demarcation threshold of positive expression and negative expression of each protein in the mutually exclusive protein pair according to the data in the negative expression cell population and the positive expression cell population comprises:

[0097] At least one of ROC and Naive Bayes methods is used to calculate the threshold value of positive expression and negative expression of each protein in the mutually exclusive protein pair; the ROC method uses the label of the positive and negative groups of the protein generated in the step c) and the expression value of the protein to generate the ROC curve and its object, and uses the Youden index to seek the optimal threshold value of the ROC curve, that is, the threshold value of positive expression and negative expression of the protein. The Naive Bayes method uses the label of the positive and negative groups of the protein generated in the step c) as the dependent variable of the model, and the expression value of the protein as the independent variable of the model to construct a Naive Bayes classification model. Then, the model is used to predict the expression value of the protein to obtain the label of the newly predicted positive and negative cells, and then the minimum value of the protein expression of the new positive group and the maximum value of the protein expression of the negative group are calculated, and the average of the minimum and maximum values ​​is taken as the threshold value of the positive and negative expression of the protein.

[0098] The COSG algorithm (Cosine Similarity-based Gene Identification) is a marker gene identification method based on cosine similarity. It is mainly used in single-cell data analysis to identify cell marker genes more accurately and quickly.

[0099] The COSG algorithm uses cosine similarity to measure the relationship between gene expression vectors. Cosine similarity evaluates the similarity between two vectors by calculating the cosine value of the angle between them in vector space. Unlike traditional statistical tests, cosine similarity compares the direction of two vectors rather than their absolute values, which means that even if two genes have the same expression pattern but different expression abundance, cosine similarity will consider them equivalent.

[0100] The COSG algorithm is particularly suitable for cell type annotation in single-cell data analysis. In single-cell sequencing technology, cell type annotation is a key step. Traditional statistical methods may encounter challenges in identifying cell marker genes because statistical tests tend to identify candidate genes with systematic differences between two groups rather than true cell markers. The COSG algorithm can more accurately identify cell marker genes by calculating the cosine similarity of gene expression vectors, thereby improving the accuracy of cell type annotation.

[0101] Compared with the existing MLTiplet CITE-seq method, this method is more stable in identifying positive and negative cell population modules. MLTiplet CITE-seq uses the mean to distinguish positive and negative cell populations, which is easily affected by background noise, clustering and expression levels. In addition, each mutually exclusive protein of MLTiplet CITE-seq needs to contain a matching gene to identify dual cell populations. This method does not require matching proteins and genes, but only requires ADT data for dual cell identification and analysis. In addition, it also provides support for background correction strategies and a variety of statistical discrimination methods, which can better adapt to different user groups.

[0102] 2-2) Identification and removal of low-quality cells: Identification and removal of low-quality cells can reduce interference with downstream data analysis. The method of the present invention integrates four low-quality cell removal strategies, including fixed threshold filtering of common quality control indicators and three public data-driven low-quality strategies (MAD, miQC and ddqc). For the fixed threshold filtering strategy, considering the variability of these indicators in different species and tissues, a fixed filtering threshold reference range is provided for each quality control parameter. These threshold ranges come from 240 single-cell studies covering 11 human tissues collected in this process, including cell UMI abundance, cell detection gene number, and cell mitochondrial gene count percentage. The threshold mode, minimum value and maximum value commonly used in these studies are further counted to provide users with a reference.

[0103] It should be pointed out that in the method of the present invention, in order to minimize the false positives of low-quality cells, this process provides a low-quality cell re-clustering process, combined with functional analysis and differential analysis, to assist in confirming cell types, so as to retain functional or rare cell types for downstream analysis. We use fixed thresholds nFeature_RNA = 200 and percent.mt = 15 to determine low-quality cells, and use this as input data for re-clustering. It can be seen that ScType and marker genes in the low-quality cell population identify the characteristics of some classic cells, such as platelets. If the user's research is on platelets, then these cells should be re-included in the analysis.

[0104] Step 3) Feature-level quality control: Features are gene or protein data in single-cell data, and each gene and protein is a feature. Feature-level quality control calculates the average expression level, expression standardization variance, and detection rate in all cells of each feature (gene / protein) in the single-cell data to be processed; corrects system background noise. Specifically includes:

[0105] 3-1) Feature (gene / protein) quality assessment: Provide a comprehensive overview of features and filter potentially low-quality gene / protein data. First, calculate the average expression level, expression standardization variance, and feature detection rate of all cells for each feature. Use visualization icons to visualize the distribution of these quality control indicators to provide a comprehensive overview of the features. To help filter potential low-quality features, this workflow provides reference thresholds for quality control indicators (the number of feature detections of all cells for each feature). These threshold ranges come from 240 single-cell studies covering 11 human tissues collected by this workflow. And further count the mode, minimum, and maximum values ​​of these thresholds to provide users with a reference. In addition, this workflow also supports the use of the scater package to calculate the degree to which each feature is affected by batch factors, assisting users in excluding features that are greatly affected by batch factors. In addition, for data from RNA+ADT combined sequencing technology, this workflow additionally provides an assessment of the correlation between RNA and corresponding surface protein features (e.g., CD3 RNA and CD3 protein).

[0106] 3-2) System background noise correction: This step can handle contamination from sources such as environmental RNA or nonspecific antibody binding to avoid affecting downstream analysis. This step integrates the DecontX and DecontPro methods to achieve global background noise correction for RNA and ADT data. For more targeted correction, this step also integrates the scCDC package to identify and remove highly contaminated RNA or protein, minimizing the risk of overcorrection of unaffected features.

[0107] Step 4) Batch effect evaluation, batch effect evaluation includes the following steps:

[0108] 4-1) Use high-dimensional dimensionality reduction technology to evaluate batch processing effects: Perform UMAP analysis on raw single-cell data at the cellular level to visualize the impact of batch effects, which is the most commonly used batch evaluation method at this stage. In addition, this process further calculates pseudo-transcriptome data for PCA analysis and aggregates count data to quickly identify batch effects between samples. This process integrates the two functions to facilitate user use, reduce the complexity of existing processes, and support more user groups.

[0109] 4-2) Evaluation of the impact of batch covariates: At this stage, the evaluation of many batch covariates is usually omitted when doing batch evaluation, and existing tools cannot analyze them well. Therefore, this process integrates singular value decomposition (SVD) analysis and the scater package to comprehensively analyze the impact of batch-related variables. Specifically, it applies SVD analysis to the pseudo-transcriptome data of single-cell data of all cells or specific cell populations at the sample level, aggregates count data to quickly identify batch-related variables that are significantly correlated with the principal components, and further checks the degree of influence of the variables to assist in judging the importance of batch-related variables. In addition, the existing publicly available scater package is used to calculate the percentage of variance explained by each batch variable for the RNA and surface protein features of the original single-cell data at the cellular level, and visualizes the contribution of these batch variables through density.

[0110] Step 5) Quality Control Web Report

[0111] In order to facilitate users to execute this process and better understand their own data sets, this process uses Rmarkdown technology to provide an integrated quality control interactive web report for the above four processes, which contains detailed parameter explanations and result descriptions, and provides a variety of interactive charts to facilitate users to quickly explore data.

[0112] The present invention also proposes a system for single cell data quality control processing based on the above method, comprising:

[0113] A sample-level quality control module is used to calculate quality control indicators in the single-cell data to be processed; identify outliers in the single-cell data to be processed based on the quality control indicators and / or cell type composition, and obtain sample-level quality control analysis results;

[0114] The cell-level quality control module is used to identify and remove empty droplets and double cell data in the single-cell data to be processed; then identify and remove low-quality cell data in the single-cell data to be processed to obtain cell-level quality control analysis results;

[0115] The feature-level quality control module is used to calculate and evaluate the average expression level, expression normalization variance and detection rate of each gene / protein in the single-cell data to be processed, and to evaluate and remove the background noise of the single-cell data to obtain the feature-level quality control analysis results;

[0116] The batch effect analysis module is used to evaluate the batch effect of the single-cell data to be processed and obtain the batch effect evaluation analysis results.

[0117] Test example

[0118] This test case performs quality control on healthy human PBMC single-cell data. The samples of healthy human PBMC single-cell data are obtained from the same healthy donor, with a sampling interval of three days, and three peripheral blood mononuclear cell (PBMCs) samples are obtained, namely TP1, TP2, and TP3. The third node sample (TP3) initially failed to pass QC, and then the data TP3-rep was remeasured. In addition, two sets of PBMC public data published on the 10x Genomics Chromium official website, namely pbmc5K and pbmc10K, were also used in the analysis. All samples used the 10x Genomics Chromium 5' gene expression platform to generate single-cell data, including RNA-seq, surface protein profiling by antibody-derived tags (ADTs), and TCR / BCR sequencing. The multi-omics fastq files generated by sequencing were processed using the 10xGenomics Cell Ranger software (v7.1.0) with the human reference genome sequence GRCh38. Among them, 4 samples (TP1, TP2, TP3, TP3-rep), with a total of 28498 cells, were included in the subsequent quality control process analysis of the method of Example 1 of the present invention.

[0119] Application Results

[0120] Some results of sample-level quality control:

[0121] First, this process can be used to roughly explore the basic situation of the multi-omics data set, such as the number of cells, the number of detected genes, the proportion of TCR chain cells detected, etc. ( Figure 2 ). In case 1, it can be found that the number of cells in the TP3 sample is relatively small and the BCR% is greater than the TCR%, which preliminarily indicates that the TP3 sample is abnormal.

[0122] For the sample abnormality module, in this case, we can find that the outlier trend of multiple quality control indicators of TP3 is found, and the outlier trend is also found in the identification of cell type ratio ( Figure 3, Table 3), furthermore, it can be concluded that the TP3 sample is abnormal. In the process of ratio comparison, the ratio of B cells and T cells in the TP3 sample deviates from the range of the reference database. Such samples are best eliminated to avoid affecting the biological authenticity.

[0123] Table 3: Cell ratios of different samples

[0124]

[0125] The identity recognition module shows that all four samples are male, which is consistent with the gender we registered and there is no error in the gender recording.

[0126] Results of cell-level quality control:

[0127] All samples can be successfully run with dual cell and low quality cell quality control methods. In this experiment, the dual cell prediction results showed that the dual cell detection ratios of different methods and different samples were different ( Figure 4-5 ).

[0128] In the determined low-quality cell re-clustering class, some cells can be seen expressing platelet characteristics such as PPBP and PF4. Such cells are over-filtered and can be further retained ( Figure 6 ).

[0129] Feature quality control part results:

[0130] At the RNA level, it was observed that cell type genes were clearly expressed among the groups. The highly contaminated genes predicted by the scCDC algorithm were not typical genes, and there was no significant improvement before and after correction, indicating that the background contamination at the RNA level was not serious and no correction was required ( Figure 7 AB). In addition, it can be seen from the ADT level that typical proteins are ubiquitously expressed in multiple cell populations, and the high contamination identified by typical proteins is greatly improved after background correction, indicating that further background contamination correction can be performed at the ADT level ( Figure 7 CD).

[0131] Some results of batch effect evaluation:

[0132] Compared with other samples, TP3-rep showed obvious batch effect in ADT sequencing, showing obvious deviation in both PCA and UMAP plots ( Figure 8 ), batch correction is required.

[0133] At the same time, the batch covariate analysis can be run smoothly. Due to the small number of samples, the p-value is not very accurate, but it can still be seen that the p-value of the antibody reagent covariate is the smallest in most pseudo-transcriptome data, and the impact on PC is very high, indicating that the antibody batch covariate has a greater impact on the protein data of most cell types than other covariates ( Fig. 9 ).

[0134] Quality control report results display:

[0135] The web quality control report can be run and displayed smoothly on a laptop. The quality control report consists of four parts. You can see the samples with problems in each chapter. The quality control report will give a yellow warning prompt and detailed explanation ( Fig.10 ).

[0136] The present invention can simultaneously perform quality control on multi-omics sequencing data of multiple samples, including scRNA-seq, surface protein sequencing, and immune repertoire sequencing, allowing users to fully evaluate their data before downstream analysis. The present invention integrates existing quality control steps, builds a bridge between various quality control processes and tools, solves the compatibility issues between different tools, and connects multiple quality control processes into one workflow. The present invention is built on the most popular Seurat framework, which is convenient for users to perform downstream analysis. In addition, the present invention is compatible with Seurat V5 and BPCells versions, can adapt to large-scale single-cell data quality control, and reduce memory usage. It has low requirements for computer memory, and an ordinary laptop can perform quality control analysis on millions of cells.

[0137] At the sample level, the present invention provides a unique outlier sample detection module, especially a cell ratio outlier detection module, which assists users in judging the biological variation of samples by integrating thousands of single-cell data samples, and assists in sample-level quality control. At the cell level, the present invention integrates 240 single-cell publications, and provides low-quality cell filtering standards and recommended reference ranges for sequencing features for 11 tissues, allowing users to select the most appropriate QC reference threshold for their experiments. In order to avoid excessive filtering of non-low-quality cells, the present invention also provides a re-clustering strategy for low-quality cells, as well as marker gene enrichment analysis. In addition, the present invention also develops a dual-cell identification strategy for ADT data in multi-omics RNA+ADT combined sequencing technology to make up for the shortcomings of dual-cell identification at the RNA level.

[0138] In addition, the present invention provides a large number of integrated and convenient calculations and QC diagrams, which greatly simplifies the complex quality control process. In particular, the present invention can generate interactive HTML reports containing single-cell data system quality control information, allowing users to quickly check data quality before downstream analysis, allowing users who do not need in-depth programming experiments to easily perform comprehensive data quality control.

[0139] Compared with other existing process tools, the single-cell multi-omics quality control process tool of the present invention is more flexible and user-friendly, providing the most comprehensive and powerful quality control workflow for single-cell multi-omics analysis.

[0140] In order to better compare the advantages and disadvantages of the methods, the present invention also tests the effect of dual cell identification using the COSG method on the ADT data in the above-mentioned RNA+ADT combined sequencing technology. The test was conducted on a total of six sets of CITE-seq data sets. Since all six sets of CITE-seq data sets contain matching VDJ sequencing data, in order to better reflect the performance difference, we selected the common B and T cell mutually exclusive protein pairs CD19-CD3 and CD20-CD3. It can be seen that the method of the present invention can effectively divide the bimodal characteristics of proteins in different data sets and proteins, while MLTiplet CITE-seq can only perform well in some samples and proteins ( Fig.11 ). Secondly, this method also maintains good segmentation performance on the corrected data of the scCDC method ( Fig.12 ).

[0141] In addition, the present invention further analyzes that the threshold fluctuations of different proteins in the iterative COSG+ statistical method strategy at different resolutions of clustering dimensions are more stable than those of MLtiplet CITE-seq ( Fig.13 ), which avoids the influence of subjective choice of resolution.

[0142] The embodiment of the present invention also provides a computer device having the above-mentioned double cell identification device. The computer device includes: one or more processors, a memory, and an interface for connecting various components, including a high-speed interface and a low-speed interface. The various components are connected to each other using different buses for communication, and can be installed on a common mainboard or installed in other ways as needed. The processor can process instructions executed in the computer device, including instructions stored in or on the memory to display graphical information of a GUI on an external input / output device (such as a display device coupled to the interface).

[0143] In some optional embodiments, if necessary, multiple processors and / or multiple buses can be used together with multiple memories and multiple memories. Similarly, multiple computer devices can be connected, and each device provides some necessary operations (for example, as a server array, a group of blade servers, or a multi-processor system). The processor can be a central processing unit, a network processor or a combination thereof. Among them, the processor can further include a hardware chip. The above-mentioned hardware chip can be a dedicated integrated circuit, a programmable logic device or a combination thereof. The above-mentioned programmable logic device can be a complex programmable logic device, a field programmable gate array, a general array logic or any combination thereof. Among them, the memory stores instructions that can be executed by at least one processor so that the at least one processor executes the double cell recognition method shown in the above embodiment.

[0144] The memory may include a program storage area and a data storage area, wherein the program storage area may store an operating system, an application required for at least one function; the data storage area may store data created according to the use of the computer device, etc. In addition, the memory may include a high-speed random access memory, and may also include a non-transient memory, such as at least one disk storage device, a flash memory device, or other non-transient solid-state storage device. In some optional embodiments, the memory may include a memory remotely arranged relative to the processor, and these remote memories may be connected to the computer device via a network. Examples of the above-mentioned network include, but are not limited to, the Internet, an intranet, a local area network, a mobile communication network, and combinations thereof.

[0145] The memory may include volatile memory, such as random access memory; the memory may also include non-volatile memory, such as flash memory, hard disk or solid state drive; the memory may also include a combination of the above types of memory.

[0146] The computer device also includes an input device and an output device. The processor, memory, input device and output device can be connected via a bus or other means. Figure 5 The example of connecting through bus is taken in the following.

[0147] The input device can receive input digital or character information, and generate key signal input related to the user settings and function control of the computer device, such as a touch screen, a keypad, a mouse, a track pad, a touch pad, an indicator bar, one or more mouse buttons, a trackball, a joystick, etc. The output device may include a display device, an auxiliary lighting device (e.g., an LED) and a tactile feedback device (e.g., a vibration motor), etc. The above-mentioned display device includes but is not limited to a liquid crystal display, a light emitting diode, a display and a plasma display. In some optional embodiments, the display device may be a touch screen.

[0148] The computer device also includes a communication interface, which is used for the computer device to communicate with other devices or a communication network.

[0149] An embodiment of the present invention also provides a computer-readable storage medium, and the above-mentioned method according to the embodiment of the present invention can be implemented in hardware, firmware, or can be implemented as a computer code that can be recorded in a storage medium, or can be implemented by downloading through a network and originally stored in a remote storage medium or a non-temporary machine-readable storage medium and will be stored in a local storage medium, so that the method described herein can be stored in such software processing on a storage medium using a general-purpose computer, a dedicated processor, or programmable or dedicated hardware.

[0150] The storage medium may be a magnetic disk, an optical disk, a read-only storage memory, a random access memory, a flash memory, a hard disk or a solid state drive, etc.; further, the storage medium may also include a combination of the above-mentioned types of memories. It is understood that a computer, a processor, a microprocessor controller or programmable hardware includes a storage component that can store or receive software or computer code. When the software or computer code is accessed and executed by the computer, processor or hardware, the method shown in the above embodiment is implemented.

[0151] Embodiments of the present application may also provide a computer program product, including computer program instructions, which, when executed by a processor, cause the processor to perform the steps in the above method. Wherein, the computer program product may be written in any combination of one or more programming languages ​​to perform program codes for performing the operations of the disclosed embodiments, wherein the programming languages ​​include object-oriented programming languages, such as Java, C++, etc., and also include conventional procedural programming languages, such as "C" language or similar programming languages. The program code may be executed entirely on a user computing device, partially on a user device, as an independent software package, partially on a user computing device, partially on a remote computing device, or entirely on a remote computing device or server.

[0152] The above embodiments are only used to illustrate the technical solutions of the embodiments of the present invention, rather than to limit them. Although the embodiments of the present invention are described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some of the technical features therein by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for quality control processing of single cell data, characterized in that: The steps include: Step 1) calculating a quality control index based on the single cell data to be processed; then identifying outliers in the single cell data to be processed based on the quality control index and / or cell type composition to obtain a sample-level quality control analysis result; Step 2) identifying empty droplets and double cell data in the single cell data to be processed; then identifying low-quality cell data in the single cell data to be processed, and drawing a quality control chart to obtain cell-level quality control analysis results; Step 3) Calculate and evaluate the average expression level, expression normalization variance and detection rate of each feature in all cells in the single-cell data to be processed to obtain feature-level quality control analysis results; Step 4) Evaluate the batch effect of the single-cell data to be processed and obtain batch effect evaluation analysis results.

2. The method for single cell data quality control processing according to claim 1, characterized in that: The single cell data is single cell multi-omics data.

3. The method for single cell data quality control processing according to claim 2, characterized in that: The single-cell multi-omics data include: scRNA-seq data, RNA+ADT sequencing data and scTCR / BCR-seq data.

4. The method for single cell data quality control processing according to claim 1, characterized in that: The quality control indicators in step 1) include: cell number, median UMI abundance of sample cells, median gene / protein detected in sample cells, median mitochondrial gene count percentage of sample cells, TCR% and BCR%.

5. The method for single cell data quality control processing according to claim 1, characterized in that: The step 1) identifying outliers in the single cell data to be processed based on the quality control index comprises the following steps: The median absolute deviation is used to calculate the deviation of the value of each quality control indicator in the processed single-cell data from the median, and the outlier samples based on each quality control indicator are identified; Preferably, the quality control indicators include the median of UMI abundance in sample cells, the median of detected genes / proteins in sample cells, and the median of mitochondrial gene count percentage in sample cells.

6. The method for single cell data quality control processing according to claim 1, characterized in that: In the step 1), identifying outliers in the single cell data to be processed based on the cell type composition comprises the following steps: 1) comparing and analyzing the cell type ratio of each sample in the single-cell data to be processed with the cell type ratio of the tissue to which the single-cell data to be processed belongs, determining whether the cell type ratio of the sample deviates from the range of the cell types in the tissue, and identifying the deviated sample as an outlier; 2) PCA combined with confidence intervals and DBSCAN method was used to identify samples with outliers in the overall cell type ratio of the sample group; Preferably, the tissues include: blood, bone, bone marrow, thymus, intestine, kidney, liver, lung, pancreas and skin; The cell types include B cells, dendritic cells, endothelial cells, fibroblasts, granulocytes, monocytes / macrophages, natural killer cells, T cells, and other cells.

7. The method for single cell data quality control processing according to claim 1, characterized in that: The step 1) further includes performing identity authentication processing on the single cell data to be processed, and identifying samples with identity authentication errors; Preferably, the identity authentication processing includes: setting a verification gene based on the sample type of the single-cell data to be processed, and then performing identity authentication according to the detection ratio of cells expressing the verification gene in the single-cell data to be processed.

8. The method for single cell data quality control processing according to claim 1, characterized in that: In the step 2), the dropletUtils package is used to identify empty droplets; Preferably, in the step 2), at least one of a single-cell transcriptome dual-cell identification method, a V(D)J sequencing dual-cell identification method, and an ADT data dual-cell identification method is used to perform dual-cell data identification; Preferably, the single-cell transcriptome dual-cell identification method includes at least one of DoubletFinder, scDbFinder, cxds, bcds, and hybrid; Preferably, the ADT data dual-cell identification method is to perform dual-cell identification based on the ADT data in the RNA+ADT combined sequencing technology by combining the COSG method with cyclic iteration.

9. The method for single cell data quality control processing according to claim 1, characterized in that: The method then identifies low-quality cell data in the single-cell data to be processed by using at least one of fixed threshold filtering, MAD filtering, miQC filtering, and ddqc filtering.

10. The method for single cell data quality control processing according to claim 1, characterized in that: In the step 2), the identified low-quality data are further subjected to cluster analysis, and the cell types of the identified low-quality data are analyzed according to the cluster analysis results; Then, based on the cell type, determine whether to retain low-quality data in the single-cell data to be processed.

11. The method for single cell data quality control processing according to claim 1, characterized in that: The step 3) further includes performing background noise correction on the single cell data to be processed, wherein the background noise correction includes: global background noise correction or local background noise correction; Preferably, the global background noise correction method includes using DecontX and DecontPro methods to perform global background noise correction on RNA and ADT data; Preferably, the local background noise correction method comprises using the scCDC method to identify and remove highly contaminating RNA or protein.

12. The method for single cell data quality control processing according to claim 1, characterized in that: The method of evaluating the batch effect of the single-cell data to be processed comprises: 1) Perform UMAP analysis on the processed single-cell data and PCA analysis on the pseudo-transcriptome data; 2) Singular value decomposition analysis and scater package are used to analyze the influence of batch-related variables.

13. The method for single cell data quality control processing according to claim 1, characterized in that: The method for single cell data quality control processing further includes: The sample-level quality control analysis results, cell-level quality control analysis results, feature-level quality control analysis results and batch effect evaluation analysis results are combined into a visual web report.

14. A system for single cell data quality control processing, characterized in that: include: Sample-level quality control module, used to calculate quality control indicators in the single-cell data to be processed; Identify outliers in the single-cell data to be processed based on the quality control indicators and / or cell type composition, and obtain sample-level quality control analysis results; The cell-level quality control module is used to identify empty droplets and double cell data in the single-cell data to be processed; then identify low-quality cell data in the single-cell data to be processed to obtain cell-level quality control analysis results; The feature-level quality control module is used to calculate and evaluate the average expression level, expression normalization variance and detection rate of each gene / protein in the single-cell data to be processed, and to evaluate and remove the background noise of the single-cell data to obtain the feature-level quality control analysis results; The batch effect analysis module is used to evaluate the batch effect of the single-cell data to be processed and obtain the batch effect evaluation analysis results.

15. An electronic device, characterized in that: comprising a memory and a processor, wherein the memory and the processor are connected; The memory stores computer instructions, and the processor executes the method for single-cell data quality control processing according to any one of claims 1 to 13 by executing the computer instructions.

16. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a computer to execute the method for single-cell data quality control processing according to any one of claims 1 to 13.

Images