Single cell transcriptome sequencing cell type automatic annotation method and device

Abstract

The invention discloses a single-cell transcriptome sequencing cell type automatic annotation method and device, relates to the technical field of single-cell transcriptome sequencing, and develops a marker gene recognition statistical method based on a built whole-organ single-cell transcriptome sequencing database and a cell marker gene bank and a Firth less-polarization logistic regression model. The marker gene score is calculated according to the preset rule and the hypothesis principle of the marker gene bank, the potential cell type of the cell cluster is identified, and the problems that an existing cell cluster marker gene identification method is high in false positive, poor in specificity, incomplete in cell marker gene bank, inaccurate in cell type annotation and the like are solved.

Description

Technical Field

[0001] The present application relates to the technical field of single-cell transcriptome sequencing, and in particular to a method and device for automatic annotation of cell types in single-cell transcriptome sequencing. Background Art

[0002] With the rapid development and widespread application of single-cell transcriptome sequencing technology, it is now possible to study gene expression characteristics at the level of a single cell, which has greatly promoted research progress in biomedical fields such as developmental biology, tumor research, immunology, and pathology. Single-cell sequencing technology allows in-depth exploration of cell heterogeneity, identification of new cell types, analysis of cell molecular regulatory networks, reconstruction of cell lineages, and exploration of biomarkers and disease subtypes.

[0003] Cell type annotation is a key step in the analysis of single-cell transcriptome sequencing data. Accurate cell type annotation is crucial for revealing cell heterogeneity, discovering new cell types, exploring cell functions and dynamic changes, and building cell maps.

[0004] Common cell type identification requires manual annotation of cell types based on genes differentially expressed in each cell cluster and searching existing literature. However, traditional differential gene identification often exhibits problems such as type I error overflow and high false positives when the proportion of cell clusters is unbalanced or gene expression data is sparse, resulting in low statistical effect. On the other hand, manual search for cell type marker genes is time-consuming and labor-intensive, and is prone to bias. The existing marker gene library is not rich in cell types and the marker genes are not accurate, which also affects the accuracy and robustness of the annotation results. Summary of the invention

[0005] In an embodiment of the present application, a method for automatic annotation of cell types in single-cell transcriptome sequencing is provided to solve the problems of low statistical effect of traditional marker gene-based annotation methods when faced with extremely unbalanced characteristics of single-cell transcriptome data, and incomplete and imperfect existing human tissue and organ cell type marker gene libraries.

[0006] In a first aspect, an embodiment of the present application provides a method for automatic annotation of cell types in single-cell transcriptome sequencing, the method comprising: inputting single-cell transcriptome sequencing data to be annotated; wherein the single-cell transcriptome sequencing data comprises a gene expression matrix and a preset cell clustering label, as well as other sample information; removing batch effects using a linear regression model; wherein the batch effect is a systematic deviation in gene expression data caused by batch factors; assigning a binary result variable value according to the cell cluster to which each cell belongs; wherein the binary result variable value indicates whether the cell belongs to a target cluster or a non-target cluster; identifying differentially expressed genes using a Firth less-biased logistic regression model; traversing any pair of target clusters and any non-target clusters, and identifying differentially expressed genes one-to-one; wherein the non-target cluster is any cell cluster other than the target cluster; calculating the average statistic of the hypothesis test of each target cluster and all other non-target clusters according to the Cauchy merging criterion; calculating the gene score according to the preset rules and the assumption principle of the single-cell marker gene library; by comparing the genes expressed in the cell clusters with the genes in the marker gene library, calculating the sum of the gene scores of overlapping genes related to a given cell type, thereby evaluating the possibility that the cell cluster belongs to a given cell type.

[0007] In a possible implementation, before inputting the single-cell transcriptome sequencing data to be annotated, the method includes: inputting a gene expression matrix and preset cell clustering labels into a marker gene library; using the gene expression matrix to identify cell clusters, evaluating the similarity between cell clusters, and merging cell clusters with high similarity.

[0008] In a possible implementation, the use of a linear regression model to remove the batch effect includes: using the batch factor as an independent variable and the gene expression level as a dependent variable in the linear regression model; fitting the linear regression model to obtain the expected expression level of each gene under the batch factor; subtracting the expected expression level from the actual expression level to obtain a residual term, and retaining the residual term as the pure gene expression level after removing the batch effect.

[0009] In a possible implementation, the method of using the Firth less-partial logistic regression model to identify differentially expressed genes includes: using a penalized likelihood function to perform parameter estimation; the expression of the penalized likelihood function is: L * (θ)=L(θ)×|I n (θ)| 0.5 Among them, L * (θ) is the penalty likelihood function, L(θ) is the standard likelihood function, which indicates the possibility of generating parameter θ from the current observation data. n (θ) is the value of the information matrix at parameter θ, |I n (θ)| 0.5is a penalty term, which represents the square root of the determinant of the information matrix. After obtaining the parameter estimation, the expression difference of each gene in different cell types is evaluated. By comparing the gene expression levels in different cell types and applying statistical tests, differentially expressed genes are identified.

[0010] In a possible implementation, the traversal of any target cluster and any non-target cluster pair, and one-to-one identification of differentially expressed genes, includes: for each target cluster, traversal of all genes, identification of the differential expression levels of genes between the target cluster and each remaining non-target cluster and completion of hypothesis testing; and obtaining a type I error of the hypothesis test for each paired target cluster and non-target group.

[0011] In a possible implementation, it also includes calculating the average statistic of multiple hypothesis tests and correcting the type I errors of multiple tests, specifically including: using the Cauchy probability value combination rule to convert multiple probability values ​​into the Cauchy statistic, and calculating the average statistic; calculating the combined type I error corresponding to the statistic according to the Cauchy distribution; and correcting the combined type I error of all genes using the Benjamin-Hochberg program.

[0012] In a possible implementation, the gene score is calculated according to preset rules and the assumption principle of the single-cell marker gene library, including: the assumption principle of the single-cell marker gene library is: excellent cell type marker genes show significant expression differences between different cell types, and the expression level of the marker gene is significantly upregulated in the target cell type; the preset rules are: differentially expressed genes are genes whose corrected type I error is less than a preset threshold; significantly upregulated differentially expressed genes are differentially expressed genes whose expression level in the target cluster is significantly higher than that in most non-target clusters; potential target cluster marker genes are significantly upregulated differentially expressed genes of the target cluster.

[0013] In a possible implementation, the following statistics matrix is ​​established for each target cluster: C G′×(M-1) Among them, G ′ is the number of significantly upregulated differentially expressed genes, M is the total number of cell clusters, M-1 is the number of non-target clusters other than the target cluster, and the elements in the matrix are the Cauchy statistics of each significantly upregulated differentially expressed gene in the differential identification between the target cluster and the non-target cluster; the gene score is calculated according to the preset rules and the assumption principle of the marker gene library, and each gene score S g for: Among them, M is the total number of cell clusters, M-1 is the number of other non-target clusters except the target cluster, and m * is the target cluster, m≠m * Indicates traversing all non-target clusters, C gm is the Cauchy statistic of gene g in cell cluster m, which is used to quantify the gene g in target cluster m* The significant difference in expression between the non-target cluster m, P gm* Indicates that in the target cluster m * The proportion of cells expressing gene g, that is, the proportion of all cells with non-zero values, P gm Represents the proportion of cells expressing gene g in the non-cell cluster m.

[0014] In a possible implementation, the method compares the genes expressed in the cell cluster with the genes in the marker gene library, calculates the sum of gene scores of overlapping genes associated with a given cell type, and thus evaluates the possibility that the cell cluster belongs to a given cell type, including: determining the genes expressed in the cell cluster, and listing a gene list according to the expressed genes; cross-comparing the gene list with the marker gene library to find overlapping genes expressed in the cell cluster and associated with a given cell type in the marker gene library; summing the gene scores of the overlapping genes to obtain a comprehensive score; wherein the comprehensive score reflects the degree of similarity between the cell cluster and the given cell type; and evaluating the possibility that the cell cluster belongs to the given cell type based on the comprehensive score; wherein the higher the comprehensive score, the higher the possibility that the cell cluster belongs to the given cell type as the comprehensive score increases.

[0015] In a second aspect, an embodiment of the present application provides an automatic annotation device for cell types in single-cell transcriptome sequencing, the device comprising: an input module for inputting single-cell transcriptome sequencing data to be annotated; wherein the single-cell transcriptome sequencing data comprises a gene expression matrix and preset cell clustering labels, as well as other sample information; a removal module for removing batch effects using a linear regression model; wherein the batch effect is a systematic deviation in gene expression data due to batch factors; an assignment module for assigning a binary result variable value according to the cell cluster to which each cell belongs; wherein the binary result variable value indicates whether the cell belongs to a target cluster or a non-target cluster; an identification module for using Firt h-less partial logistic regression model to identify differentially expressed genes; a traversal module, used to traverse any target cluster and any non-target cluster pair, and identify differentially expressed genes one-to-one; wherein the non-target cluster is any cell cluster other than the target cluster; a first calculation module, used to calculate the average statistic of the hypothesis test of each target cluster and all other non-target clusters according to the Cauchy merging criterion; a second calculation module, used to calculate the gene score according to preset rules and the assumption principle of the single-cell marker gene library; an evaluation module, used to calculate the sum of the gene scores of overlapping genes related to a given cell type by comparing the genes expressed in the cell cluster with the genes in the marker gene library, thereby evaluating the possibility that the cell cluster belongs to a given cell type.

[0016] One or more technical solutions provided in the embodiments of the present application have at least the following technical effects: the embodiments of the present application provide a method for automatic annotation of cell types for single-cell transcriptome sequencing, which uses a linear regression model to remove batch effects, and can eliminate the systematic bias that may be introduced during the processing of different batches of samples, thereby more accurately reflecting the gene expression differences between cells. The Firth less partial logistic regression model is used for differential gene detection, which can more accurately identify genes that are differentially expressed between different cell clusters, thereby improving the accuracy of cell type annotation. The possibility that a cell cluster belongs to a specific cell type is evaluated by calculating the sum of gene scores of overlapping genes associated with a given cell type. This method not only takes into account the expression differences of a single gene, but also comprehensively considers the combined effects of multiple genes, thereby enhancing the reliability and persuasiveness of data interpretation. This method is not only applicable to single-cell transcriptome sequencing data, but can also be extended to other types of omics data, such as single-cell epigenetic data, single-cell proteomic data, etc. Therefore, it has wide applicability and application prospects. It solves the problem that traditional annotation methods based on marker genes often encounter problems when faced with extremely unbalanced characteristics of single-cell transcriptome data, such as model parameter hypothesis test statistics deviating from normal distribution and type I errors (α errors) being difficult to control, resulting in deviations in cell type judgments and affecting the accuracy of annotation results. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0018] Figure 1 A flowchart of a method for automatic annotation of cell types in single-cell transcriptome sequencing provided in an embodiment of the present application;

[0019] Figure 2 A specific flow chart of using a linear regression model to remove batch effects provided in an embodiment of the present application;

[0020] Figure 3 A specific flow chart for identifying differentially expressed genes using the Firth less-biased logistic regression model provided in the embodiments of the present application;

[0021] Figure 4 Provided for the embodiments of the present application are also specific flow charts for calculating the average statistic of multiple hypothesis tests and correcting multiple test type I errors;

[0022] Figure 5A specific flow chart for calculating the sum of gene scores of overlapping genes associated with a given cell type, thereby evaluating the possibility that a cell cluster belongs to a given cell type, provided in an embodiment of the present application;

[0023] Figure 6 A schematic diagram of an automatic annotation device for single-cell transcriptome sequencing cell types provided in an embodiment of the present application;

[0024] Figure 7 Schematic diagram of the single-cell transcriptome sequencing cell type automatic annotation server provided in the embodiments of the present application. DETAILED DESCRIPTION

[0025] The technical solutions in the embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0026] The following describes some of the techniques involved in the embodiments of the present application to facilitate understanding, and they should be considered as merely exemplary. Therefore, it should be appreciated by those of ordinary skill in the art that various changes and modifications may be made to the embodiments described herein without departing from the scope and spirit of the present application. Similarly, for the sake of clarity and conciseness, some descriptions of well-known functions and structures are omitted in the following description.

[0027] The present application embodiment provides a method for automatic annotation of cell types in single-cell transcriptome sequencing, such as Figure 1 As shown, the method includes steps S101 to S108. Among them, Figure 1 This is only an execution order shown in the embodiment of the present application, and does not represent the only execution order of a method for automatic annotation of single-cell transcriptome sequencing cell types. If the final result can be achieved, Figure 1 The steps shown may be performed in parallel or reversed.

[0028] S101: Input the single-cell transcriptome sequencing data to be annotated, wherein the single-cell transcriptome sequencing data includes a gene expression matrix and preset cell clustering labels, as well as other sample information.

[0029] Specifically, before inputting the single-cell transcriptome sequencing data to be annotated, the following steps are included: inputting a gene expression matrix and preset cell clustering labels into a labeling gene library (CellLabeler). Using the gene expression matrix, cell clusters are identified, similarities between cell clusters are evaluated, and cell clusters with high similarity are merged.

[0030] Specifically, the gene expression matrix in the present application is a logarithmically normalized gene expression matrix, and its expression is: G×N . Where G represents the number of genes and N represents the number of cells. Each row of the gene expression matrix represents a gene, with a total of G genes, and each column represents a cell, with a total of N cells. The gene expression matrix can be regarded as a table, and the value in each cell represents the expression level of the corresponding gene in the corresponding cell. Some prior information or preliminary clustering results can be input into the marker gene library as pre-clustering labels. The expression of the pre-clustering label is: L N = {l i |l i =1,2,…,M; i =1,2,…,N}. Among them, L N is a set, representing the pre-clustering labels of all cells. The subscript N indicates that this set contains N elements, that is, there are N cells in the data set. i For the set L N An element in represents the pre-clustering label of the ith cell. M is the number of predefined cell clusters. i The value range is an integer from 1 to M, indicating that the i-th cell is classified as one of the M possible cell clusters. i is an index variable used to identify each cell in the data set. The pre-clustering label can tell which cell cluster each cell is classified into during the clustering process. Based on the gene expression matrix, unsupervised learning methods (such as hierarchical clustering or clustering algorithms) are used to identify different cell clusters. These cell clusters represent cell populations that are similar at the transcriptome level. In the clustering process, it is usually necessary to adjust some parameters (such as the number of clusters, similarity threshold, etc.) according to the specific situation of the data and the experimental goals to obtain the best results.

[0031] Specifically, after identifying the cell clusters, the similarities between these cell clusters can be evaluated by calculating statistics between different cell clusters (such as correlation of gene expression, distance metric, etc.). Dimensionality reduction techniques (such as PCA) are used to project high-dimensional gene expression data into two-dimensional or three-dimensional space in order to more intuitively evaluate the similarities between cell clusters.

[0032] Furthermore, in single-cell data analysis, it is a common strategy to compare the cell cluster of interest (target cluster) with all other cell clusters one by one, a strategy called "pairwise" comparison. When only the cell cluster of interest is compared with the "other" or "remaining" cells as a whole, the analysis results will be distorted due to data imbalance (i.e., the number of cells in the target cluster is very different from the number of cells in all other clusters). Therefore, it is necessary to set a threshold to determine which cell clusters should be merged based on the results of the similarity assessment. This threshold needs to be determined based on the specific circumstances of the data and the experimental goals. For cell clusters with a similarity higher than the threshold, they can be merged into a larger cluster. This can be done by rerunning the clustering algorithm or using other methods.

[0033] S102: Using a linear regression model to remove batch effects, wherein batch effects are systematic deviations in gene expression data caused by batch factors.

[0034] Specifically, batch factors refer to the fact that samples are divided into different batches or groups due to differences in experimental conditions, time, reagent batches, etc. in experiments or data analysis. Other factors include changes in environmental factors such as temperature, humidity, and light, as well as sample collection time, storage conditions, and processing methods.

[0035] S103: assigning a binary result variable value according to the cell cluster to which each cell belongs, wherein the binary result variable value indicates whether the cell belongs to the target cluster or the non-target cluster.

[0036] Figure 2 A specific flow chart of using a linear regression model to remove batch effects is provided in the embodiment of the present application, such as Figure 2 As shown, it includes steps S201 to S203.

[0037] S201: In the linear regression model, batch factors were used as independent variables and gene expression levels were used as dependent variables.

[0038] Specifically, in the linear regression model, the batch factor is used as the independent variable to quantify the degree of influence of the batch effect on gene expression. The model can predict the expected value of gene expression under a given batch factor. Gene expression refers to the number of mRNA molecules transcribed from a gene in a cell. In the linear regression model, gene expression is the dependent variable and is the target to be explained and predicted. Through the model, we can understand the impact of batch factors on gene expression and try to eliminate this impact to obtain a purer biological difference.

[0039] S202: Fitting a linear regression model to obtain the expected expression level of each gene under the batch factor.

[0040] Specifically, in this step, a linear regression model is fitted using the collected batch factors and the corresponding gene expression data. The goal of this model is to describe how batch factors affect gene expression. Through least squares or other optimization methods, the parameters of the model (such as slope, intercept, etc.) can be obtained, which define the linear relationship between batch factors and gene expression. Once the model is fitted, the model can be used to predict the expected expression of each gene under a given batch factor. These expected expression levels are calculated based on batch factors rather than biological differences, so they represent the impact of batch effects on gene expression.

[0041] S203: Subtract the expected expression amount from the actual expression amount to obtain a residual term, and retain the residual term as the pure gene expression amount after removing the batch effect.

[0042] Specifically, in this step, the actual expression of each gene (i.e., the original measurement value) is subtracted from the expected expression predicted by the model to obtain the residual term. The residual term represents the difference between the actual expression and the expected expression, which is due to biological differences rather than batch effects. Therefore, the residual term can be regarded as the pure gene expression after removing the batch effect. These residual terms retain the biological difference information while excluding the influence of batch effects.

[0043] Specifically, it is necessary to determine which cell cluster is of interest (target cluster) based on the research question, and which cell cluster or clusters will be used as reference (non-target cluster). For each cell, a binary outcome variable value is assigned according to the cell cluster to which it belongs. 1 can be used to represent cells belonging to the target cluster, and 0 can be used to represent cells belonging to the non-target cluster.

[0044] S104: Firth less-partial logistic regression model was used to identify differentially expressed genes.

[0045] Specifically, as the imbalance of gene expression (dependent variable) (such as sparse sample distribution, quasi-complete separation, zero overfilling, extreme class imbalance, etc.) becomes stronger, the bias of prediction probability will become more serious, so the Firth biased logistic regression model is used for differential gene detection. Among them, the balanced correction logistic regression model is Firth's bias-corrected logistic regression model.

[0046] Figure 3 The specific flow chart of using Firth's less partial logistic regression model to identify differentially expressed genes provided in the embodiment of the present application is as follows: Figure 3 As shown, it includes steps S301 to S303.

[0047] S301: Use the penalized likelihood function to perform parameter estimation. The expression of the penalized likelihood function is: L *(θ)=L(θ)×|I n (θ)| 0.5 Among them, L * (θ) is the penalty likelihood function, L(θ) is the standard likelihood function, which indicates the possibility of generating parameter θ from the current observation data. n (θ) is the value of the information matrix at parameter θ, |I n (θ)| 0.5 is the penalty term, which represents the square root of the determinant of the information matrix.

[0048] Specifically, the information matrix is ​​a second-order derivative matrix that contains information about the uncertainty of the parameter estimates. n (θ)| 0.5 is the square root of the information matrix determinant, which is introduced into the penalized likelihood function as a penalty term. The purpose of this penalty term is to adjust the parameter estimate to make it more stable. Specifically, when the determinant of the information matrix is ​​small, that is, when the uncertainty of the parameter estimate is large, the penalty term will increase, thereby reducing the amplitude of the parameter estimate and making it more stable.

[0049] S302: After obtaining parameter estimates, evaluate the expression differences of each gene in different cell types.

[0050] Specifically, based on the obtained parameter estimates, the balanced corrected logistic regression model can now predict the expression level of each gene in different cell types. These predicted values ​​will be used to compare the expression of genes in different cell types.

[0051] S303: Compare gene expression levels in different cell types and apply statistical tests to identify differentially expressed genes.

[0052] Specifically, differentially expressed genes refer to genes whose expression levels differ significantly between different cell types. Commonly used statistical test methods include analysis of variance, etc. The test method calculates a statistic, and then based on the statistic and a pre-set significance level, determines whether the expression difference of the gene in different cell types is significant.

[0053] S105: Traverse any pair of target clusters and any non-target clusters to identify differentially expressed genes one by one, wherein the non-target cluster is any cell cluster other than the target cluster.

[0054] Traverse any target cluster and any non-target cluster pair, and identify differentially expressed genes one by one, including: for each target cluster, traverse all genes, identify the differential expression level of genes between the target cluster and each remaining non-target cluster and complete hypothesis testing. Obtain the type I error of the hypothesis test for each paired target cluster and non-target group.

[0055] Specifically, a difference analysis is performed on each gene between the target cluster and the non-target cluster to obtain multiple probability values, including: for each target cluster, a specific gene is selected and its expression data in the target cluster is compared with the expression data of this gene in other non-target clusters. The comparison is completed by applying statistical tests, and each statistical test will produce a probability value. Among them, the probability value is the probability that the observed difference is caused by random errors alone.

[0056] Furthermore, for a target cluster m and a gene g, in the process of differential analysis, M-1 p values, i.e., probability values, of the gene g in the target cluster m and M-1 non-target cell clusters will be obtained. Each p value represents the probability that the observed difference is only caused by random error under the null hypothesis (i.e., there is no difference in the expression of gene g between the target cluster m and the corresponding non-target cluster).

[0057] Specifically, it is first necessary to determine the target cluster and one or more non-target clusters of the study. For each target cluster, select a specific gene from all the genes studied. Collect the expression data of the gene in the target cluster, and also collect the expression data of the gene in all non-target clusters. Compare the gene expression data in the target cluster with the gene expression data in the non-target cluster. This step usually involves statistical tests to quantify the differences between the two sets of data. A probability value (p-value) can be obtained through statistical tests. This probability value indicates the possibility that the observed differences are caused only by random errors. If the p-value is very small (usually a threshold is set, such as p<0.05), it is believed that the observed differences are not caused by random errors, but are real.

[0058] S106: Calculate the average statistic of the hypothesis test between each target cluster and all other non-target clusters according to the Cauchy merging criterion.

[0059] Figure 4 The specific flow chart provided for the embodiment of the present application also includes calculating the average statistic of multiple hypothesis tests and correcting multiple test type I errors, such as Figure 4 As shown, it includes steps S401 to S403.

[0060] S401: convert multiple probability values ​​into Cauchy statistics using the Cauchy probability value combination rule, and calculate the average statistic.

[0061] Specifically, multiple p-values ​​(probability values) obtained by comparing all non-target clusters (reference clusters) with each gene in the target cluster are combined. For example, if there are 3 non-target clusters and 1 target cluster, then for each gene in the target cluster, 3 p-values ​​will be obtained, and these p-values ​​need to be combined.

[0062] S402: Calculate the combined type I error corresponding to the statistic according to the Cauchy distribution.

[0063] S403: The combined type I error of all genes was corrected using the Benjamin-Hochberg procedure.

[0064] Specifically, the Cauchy probability value combination rule is the Cauchy P value combination rule. The Cauchy P value combination rule is used to integrate these M-1 p values. First, each p value obtained by the balanced correction logistic regression model is converted into a corresponding Cauchy statistic. The Cauchy statistic is a function based on the p value, and its characteristics can better handle the extreme value problem that may occur when the p value is combined. Then, these converted Cauchy statistics are summed. The summation here does not directly add the probability values, but first converts each probability value into a Cauchy statistic. The Cauchy statistic is a function based on the probability value, which is used to provide better statistical properties when combining probability values. Then, all these converted Cauchy statistics are added.

[0065] Specifically, the summed Cauchy statistic is converted back to a single probability value, i.e., a single p-value, based on the standard Cauchy distribution, which is the same standard probability value as the number of genes. This single p-value combines the information of the original M-1 p-values ​​and provides a comprehensive measure of the expression difference of gene g between the target cluster m and the non-target cluster.

[0066] Specifically, the summation result is a statistic that combines multiple original probability values. This step converts the summation result back into a standard probability value, which is calculated based on the standard Cauchy distribution. Due to the characteristics of the Cauchy distribution, the standard probability value obtained in this step can better reflect the combined effect of the original multiple probability values.

[0067] Specifically, the last step is to correct the standard probability value to control the false positive rate caused by multiple statistical tests. The Bayesian-Yekutieli procedure is a method for controlling the false positive discovery rate (FDR). The procedure provides a corrected standard probability value for each gene, which can more accurately reflect the true significance of the difference in expression between the target cluster m and the non-target cluster.

[0068] S107: Calculate gene scores based on preset rules and assumptions of the single-cell marker gene library.

[0069] The assumption principle of the single-cell marker gene library is that excellent cell type marker genes show significant expression differences between different cell types, and the expression level of marker genes is significantly upregulated in the target cell type.

[0070] Specifically, good cell type marker genes should show significant expression differences between different cell types. This means that these genes are highly expressed in the cell type of interest (target cell type) and low or no expression in other cell types. Upregulated expression is a key feature of cell type specificity and helps determine which genes can serve as markers for specific cell types.

[0071] The preset rule is: differentially expressed genes are genes whose corrected type I errors are less than the preset threshold. Significantly up-regulated differentially expressed genes are differentially expressed genes whose expression levels in the target cluster are significantly higher than those in most non-target clusters. Potential target cluster marker genes are significantly up-regulated differentially expressed genes in the target cluster. Among them, the preset threshold is 0.05, and most refers to more than 80% of the clusters in the non-target clusters, which means that it is necessary to compare the expression levels of genes in the target cell cluster with the expression levels of genes in the non-target cell clusters, and ensure that the expression of the gene is low in most non-target clusters.

[0072] Genes that meet the preset rules are collected into a new gene set. The gene expression differences between each target cluster and non-target cluster are compared to form a probability value matrix. Each element in the probability value matrix represents the probability value of the expression difference between the target cluster and the non-target cluster.

[0073] Specifically, the following statistical matrix is ​​established for each target cluster: C G′×(M-1) Among them, G ′ is the number of significantly upregulated differentially expressed genes, M is the total number of cell clusters, M-1 is the number of non-target clusters other than the target cluster, and the elements in the matrix are the Cauchy statistics of each significantly upregulated differentially expressed gene in the differential identification of the target cluster and the non-target cluster. Each element in the probability value matrix C[i][j] represents the probability value of the expression difference between the i-th gene in the j+1-th non-target cluster (because the index usually starts from 0) and the target cluster.

[0074] Furthermore, calculating the gene score according to the preset rules and the assumption principle of the marker gene library also includes: calculating the gene score in two ways.

[0075] The first way to calculate the gene score is Among them, M is the total number of cell clusters, M-1 is the number of other non-target clusters except the target cluster, and m * is the target cluster, m≠m * Indicates traversing all non-target clusters, C gm is the Cauchy statistic of gene g in cell cluster m, which is used to quantify the gene g in target cluster m * The significant difference in expression between the non-target cluster m and Indicates that in the target cluster m * The proportion of cells expressing gene g, that is, the proportion of all cells with non-zero values, P gm Represents the proportion of cells expressing gene g in the non-cell cluster m.

[0076] The second way to calculate gene scores is: Among them, S g is the gene score, is the gene g in the target cluster m * The Cauchy statistic in , M is the total number of cell clusters, M-1 is the number of non-target clusters other than the target cluster, m is the cell cluster, m * is the target cluster, m≠m * Indicates that only non-target clusters are considered. Indicates that in the target cluster m * The proportion of cells with gene g in gm Represents the proportion of cells with gene g in candidate cell cluster m.

[0077] S108: By comparing the genes expressed in the cell cluster with the genes in the marker gene library, the sum of gene scores of overlapping genes associated with a given cell type is calculated to assess the possibility that the cell cluster belongs to a given cell type.

[0078] Figure 5 A specific flow chart for calculating the sum of gene scores of overlapping genes associated with a given cell type, thereby evaluating the possibility that a cell cluster belongs to a given cell type, is provided in an embodiment of the present application, such as Figure 5 As shown, it includes steps S501 to S504.

[0079] S501: Determine the genes expressed in the cell clusters and make a gene list based on the expressed genes.

[0080] Specifically, it is necessary to extract all expressed genes from the sequencing data of the cell cluster and generate a gene list.

[0081] S502: Cross-match the gene list with the marker gene library to find overlapping genes that are expressed in the cell clusters and are associated with a given cell type in the marker gene library.

[0082] Specifically, the gene list obtained in the previous step is compared with the marker gene library. The marker gene library is a database containing known cell type-specific marker genes. By comparing, overlapping genes that are expressed in the cell cluster and are associated with a given cell type in the marker gene library can be found.

[0083] S503: summing up the gene scores of the overlapping genes to obtain a comprehensive score, wherein the comprehensive score reflects the degree of similarity between the cell cluster and the given cell type.

[0084] S504: Evaluate the possibility that the cell cluster belongs to a given cell type according to the comprehensive score, wherein the higher the comprehensive score is, the higher the possibility that the cell cluster belongs to the given cell type is.

[0085] Specifically, this can be achieved by setting a judgment threshold. If the comprehensive score exceeds the judgment threshold, it is considered that the cell cluster belongs to a given cell type.

[0086] Furthermore, the marker gene library of this application is the CellLabeler database, which contains cell type-specific marker genes from a variety of common healthy human organs and tissues. These cell type annotations and corresponding marker genes are carefully collected, standardized and integrated from a large number of authoritative single-cell RNA sequencing studies. And in this application, genes with highly specific expression identified in cell types that are less mentioned in existing literature are also included in the CellLabeler database after sufficient verification and demonstration.

[0087] The present application embodiment also provides a single-cell transcriptome sequencing cell type automatic annotation device 600, such as Figure 6 As shown, the device includes: an input module 601, a removal module 602, an assignment module 603, an identification module 604, a traversal module 605, a first calculation module 606, a second calculation module 607 and an evaluation module 608.

[0088] The input module 601 is used to input the single-cell transcriptome sequencing data to be annotated, wherein the single-cell transcriptome sequencing data includes a gene expression matrix and preset cell clustering labels, as well as other sample information.

[0089] The removal module 602 is used to remove the batch effect using a linear regression model, wherein the batch effect is a systematic deviation in the gene expression data caused by batch factors.

[0090] The assignment module 603 is used to assign a binary result variable value according to the cell cluster to which each cell belongs, wherein the binary result variable value indicates whether the cell belongs to the target cluster or the non-target cluster.

[0091] The identification module 604 is used to identify differentially expressed genes using the Firth less-partial logistic regression model.

[0092] The traversal module 605 is used to traverse any target cluster and any non-target cluster pair, and identify differentially expressed genes one-to-one. The non-target cluster is any cell cluster other than the target cluster.

[0093] The first calculation module 606 is used to calculate the average statistic of the hypothesis test between each target cluster and all other non-target clusters according to the Cauchy merging criterion.

[0094] The second calculation module 607 is used to calculate the gene score according to the preset rules and the assumption principle of the single cell marker gene library.

[0095] The evaluation module 608 is used to compare the genes expressed in the cell cluster with the genes in the marker gene library, calculate the sum of gene scores of overlapping genes related to a given cell type, and thus evaluate the possibility that the cell cluster belongs to a given cell type.

[0096] Some modules in the apparatus described in the present application can be described in the general context of computer executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, classes, etc. that perform specific tasks or implement specific abstract data types. The present application can also be practiced in distributed computing environments, in which tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules can be located in local and remote computer storage media including storage devices.

[0097] The devices or modules described in the above application embodiments can be implemented by computer chips or entities, or by products with certain functions. For the convenience of description, the above devices are described in various modules according to their functions. When implementing the embodiments of the present application, the functions of each module can be implemented in the same or multiple software and / or hardware. Of course, the module that implements a certain function can also be implemented by combining multiple sub-modules or sub-units.

[0098] The methods, devices or modules described in this application can be implemented in a computer-readable program code mode. The controller can be implemented in any appropriate manner. For example, the controller can take the form of a microprocessor or processor and a computer-readable medium storing a computer-readable program code (such as software or firmware) that can be executed by the (micro) processor, a logic gate, a switch, an application-specific integrated circuit (English: Application Specific Integrated Circuit, referred to as: ASIC), a programmable logic controller and an embedded microcontroller. Examples of controllers include but are not limited to the following microcontrollers: ARC 625D, Atmel AT91SAM, Microchip PIC18F26K20 and Silicone Labs C8051F320. The memory controller can also be implemented as part of the control logic of the memory. Those skilled in the art also know that in addition to implementing the controller in a pure computer-readable program code mode, the controller can be completely implemented in the form of a logic gate, a switch, an application-specific integrated circuit, a programmable logic controller and an embedded microcontroller by logically programming the method steps. Therefore, this controller can be considered as a hardware component, and the devices included in it for implementing various functions can also be regarded as structures within the hardware component. Or even, the means for realizing various functions may be regarded as both a software module for realizing the method and a structure within a hardware component.

[0099] like Figure 7 As shown, an embodiment of the present application also provides a single-cell transcriptome sequencing cell type automatic annotation server, including a memory 701 and a processor 702; the memory 701 is used to store computer-executable instructions; the processor 702 is used to execute computer-executable instructions to implement a single-cell transcriptome sequencing cell type automatic annotation method described above in the embodiment of the present application.

[0100] An embodiment of the present application also provides a computer-readable storage medium, which stores executable instructions. When a computer executes the executable instructions, it can implement the method for automatic annotation of cell types in single-cell transcriptome sequencing described above in the embodiment of the present application.

[0101] It can be known from the description of the above implementation mode that a person skilled in the art can clearly understand that the present application can be implemented by means of software plus necessary hardware. Based on such an understanding, the technical solution of the present application can be essentially or partly contributed to the prior art in the form of a software product, or it can be reflected in the implementation process of data migration. The computer software product can be stored in a storage medium, such as ROM / RAM, a magnetic disk, an optical disk, etc., including a number of instructions to enable a computer device (which can be a personal computer, a mobile terminal, a server, or a network device, etc.) to execute the method described in the embodiment of the present application.

[0102] The various embodiments in this specification are described in a progressive manner, and the same or similar parts between the various embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. All or part of this application can be used in many general or special computer system environments or configurations.

[0103] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit the present application. Although the present application has been described in detail with reference to the aforementioned embodiments, a person of ordinary skill in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some or all of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the present application.

Claims

1. A method for automatic annotation of cell types in single-cell transcriptome sequencing, characterized in that: include: Input the single-cell transcriptome sequencing data to be annotated; wherein the single-cell transcriptome sequencing data includes a gene expression matrix and preset cell clustering labels, as well as other sample information; A linear regression model was used to remove batch effects, where batch effects are systematic biases in gene expression data caused by batch factors. Assigning a binary result variable value according to the cell cluster to which each cell belongs; wherein the binary result variable value indicates whether the cell belongs to the target cluster or the non-target cluster; Firth less-partial logistic regression model was used to identify differentially expressed genes; Traversing any pair of target cluster and any non-target cluster, identifying differentially expressed genes one-to-one; wherein the non-target cluster is any cell cluster other than the target cluster; According to the Cauchy merging criterion, the average statistic of the hypothesis test between each target cluster and all other non-target clusters is calculated; Gene scores are calculated based on preset rules and assumptions of single-cell marker gene libraries; The likelihood that a cell cluster belongs to a given cell type is assessed by comparing the genes expressed in the cell cluster with the genes in the marker gene library and calculating the sum of gene scores of overlapping genes associated with a given cell type.

2. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: Before the single-cell transcriptome sequencing data to be annotated is input, it includes: Input the gene expression matrix and preset cell cluster labels into the marker gene library; The gene expression matrix is ​​used to identify cell clusters, evaluate the similarity between cell clusters, and merge cell clusters with high similarity.

3. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: The linear regression model is used to remove batch effects, including: In the linear regression model, batch factors were used as independent variables and gene expression levels as dependent variables; Fit the linear regression model to obtain the expected expression of each gene under batch factors; The actual expression was subtracted from the expected expression to obtain the residual term, which was retained as the pure gene expression after removing the batch effect.

4. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: The method of using the Firth less-partial logistic regression model to identify differentially expressed genes includes: Use penalized likelihood function for parameter estimation; The expression of the penalty likelihood function is: L * (θ)=L(θ)×|I n (θ)| 0.5 Among them, L * (θ) is the penalty likelihood function, L(θ) is the standard likelihood function, which indicates the possibility of generating parameter θ from the current observation data. n (θ) is the value of the information matrix at parameter θ, |I n (θ)| 0.5 is the penalty term, which represents the square root of the determinant of the information matrix; After obtaining parameter estimates, the expression differences of each gene in different cell types were evaluated; By comparing gene expression levels in different cell types, differentially expressed genes were identified using statistical tests.

5. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: The traversing any pair of target clusters and any non-target clusters to identify differentially expressed genes one by one includes: For each target cluster, all genes are traversed to identify the differential expression levels of genes between the target cluster and each of the remaining non-target clusters and complete hypothesis testing; The Type I error of the hypothesis test is obtained for each pair of target cluster and non-target group.

6. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 5, characterized in that: It also includes calculating the average statistic for multiple hypothesis tests and correcting for multiple test type I errors, including: Use the Cauchy probability value combination rule to convert multiple probability values ​​into Cauchy statistics and calculate the average statistic; The pooled type I error corresponding to this statistic is calculated according to the Cauchy distribution; The Benjamin-Hochberg procedure was used to correct for the pooled type I error for all genes.

7. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: The gene score is calculated according to the preset rules and the assumption principle of the single cell marker gene library, including: The assumption principle of the single-cell marker gene library is that excellent cell type marker genes show significant expression differences between different cell types, and the expression level of the marker gene is significantly upregulated in the target cell type; The preset rule is: differentially expressed genes are genes whose corrected type I errors are less than a preset threshold; Significantly up-regulated differentially expressed genes are differentially expressed genes whose expression levels in the target cluster are significantly higher than those in most non-target clusters; The potential target cluster marker genes are the significantly upregulated differentially expressed genes of the target cluster.

8. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 7, characterized in that: For each target cluster, the following statistical matrix is ​​established: C m′×(M-1) Among them, G ′ is the number of significantly upregulated differentially expressed genes, M is the total number of cell clusters, M-1 is the number of non-target clusters other than the target cluster, and the elements in the matrix are the Cauchy statistics of each significantly upregulated differentially expressed gene in the differential identification between the target cluster and the non-target cluster; The gene score is calculated according to the preset rules and the assumption principle of the marker gene library. Each gene score S g for: Among them, M is the total number of cell clusters, M-1 is the number of non-target clusters other than the target cluster, and m * is the target cluster, m≠m * Indicates traversing all non-target clusters, C gm is the Cauchy statistic of gene g in cell cluster m, which is used to quantify the gene g in target cluster m * The significant difference in expression between the non-target cluster m and Indicates that in the target cluster m * The proportion of cells expressing gene g in the expression, that is, the proportion of all cells with non-zero values, P gm Represents the proportion of cells expressing gene g in the non-cell cluster m.

9. The method for automatic annotation of cell types by single-cell transcriptome sequencing according to claim 1, characterized in that: The method compares the genes expressed in the cell cluster with the genes in the marker gene library, calculates the sum of gene scores of overlapping genes associated with a given cell type, and thereby evaluates the possibility that the cell cluster belongs to a given cell type, including: Determine the genes expressed in the cell clusters and make a gene list based on the genes expressed; Cross-matching the gene list with the marker gene library to find overlapping genes expressed in the cell clusters and associated with a given cell type in the marker gene library; The gene scores of the overlapping genes are summed to obtain a comprehensive score; wherein the comprehensive score reflects the degree of similarity between the cell cluster and the given cell type; The possibility that the cell cluster belongs to a given cell type is evaluated according to the comprehensive score; wherein, the higher the comprehensive score is, the higher the possibility that the cell cluster belongs to the given cell type is.

10. A single-cell transcriptome sequencing cell type automatic annotation device, characterized in that: include: An input module, used to input single-cell transcriptome sequencing data to be annotated; wherein the single-cell transcriptome sequencing data includes a gene expression matrix and preset cell clustering labels, as well as other sample information; A removal module is used to remove batch effects using a linear regression model; batch effects are systematic deviations in gene expression data caused by batch factors; An assignment module, used for assigning a binary result variable value according to the cell cluster to which each cell belongs; wherein the binary result variable value indicates whether the cell belongs to a target cluster or a non-target cluster; Identification module, used to identify differentially expressed genes using Firth's less partial logistic regression model; A traversal module, used to traverse any pair of target clusters and any non-target clusters, and identify differentially expressed genes one-to-one; wherein the non-target cluster is any cell cluster other than the target cluster; The first calculation module is used to calculate the average statistic of the hypothesis test of each target cluster and all other non-target clusters according to the Cauchy merging criterion; The second calculation module is used to calculate the gene score according to the preset rules and the assumption principle of the single-cell marker gene library; The evaluation module is used to evaluate the possibility that the cell cluster belongs to a given cell type by comparing the genes expressed in the cell cluster with the genes in the marker gene library and calculating the sum of gene scores of the overlapping genes related to the given cell type.

Images