Animal single cell data cell type annotation method and system

Abstract

The invention discloses an animal single cell data cell type annotation method, which is characterized by comprising the following steps: collecting Bulk RNA-seq data of purified cell types in various tissues and organs of a specific animal, and integrating the Bulk RNA-seq data into a reference data set after preprocessing; the method comprises the following steps: acquiring single-cell RNA-seq original data, screening the single-cell RNA-seq original data to obtain high-quality cells, screening high-variation genes from the high-quality cells, processing the high-variation genes, extracting principal components, performing dimensionality reduction on the principal components, and performing cell clustering based on a dimensionality reduction result to obtain cell clusters; and calculating expression similarity between the cell clusters and the reference data set based on the high-variation genes, determining initial cell types of the cell clusters according to a similarity result, and carrying out iterative tuning on the initial cell types with similar scores to obtain a final cell type annotation result. The method provides efficient and accurate technical support for animal single cell research.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of single-cell transcriptomics, and particularly relates to an animal single-cell data cell type annotation method and system. BACKGROUND

[0002] In recent years, the rapid development of single-cell RNA sequencing (scRNA-seq) technology has enabled researchers to analyze gene expression heterogeneity in complex tissues at single-cell resolution. This technology captures the transcriptome information of individual cells, revealing the composition, state changes and functional diversity of cell populations, providing an unprecedented research perspective for developmental biology, immunology, oncology and neuroscience. At present, mainstream single-cell sequencing platforms have achieved high-throughput and high-sensitivity single-cell transcriptome analysis, making large-scale single-cell mapping possible.

[0003] At the data analysis level, the processing flow of scRNA-seq data mainly includes quality control (QC), data standardization, dimensionality reduction (such as PCA, t-SNE, UMAP), cell clustering and differential expression analysis. Common analysis tools (such as Seurat, Scanpy, Monocle, etc.) provide a complete computational framework to help researchers extract biological information from raw sequencing data. Among them, cell clustering is one of the core steps, which divides cells into different subgroups by calculating the similarity of gene expression between cells. However, the clustering result only reflects the similarity of gene expression patterns, and cannot directly correspond to known cell types. Therefore, how to accurately annotate the biological identity of these cell clusters has become one of the key challenges of single-cell data analysis.

[0004] Although single-cell technology has made significant progress, cell type annotation still faces many problems. First, current methods are highly dependent on known marker genes for manual annotation. Researchers often need to consult literature or databases to obtain marker genes for specific cell types, and based on the expression pattern of these genes, the clustering results are manually annotated. This process not only consumes time and effort, but also easily introduces subjective bias, especially when dealing with complex tissues or rare cell types, marker genes may not be specific enough, leading to annotation errors or omissions. In addition, different research teams may use inconsistent marker gene standards, making it difficult to compare annotation results across datasets.

[0005] Secondly, existing annotation methods generally ignore species differences, and most tools and reference databases are mainly developed based on human or mouse data. Due to significant differences in gene expression profiles between different species, directly applying human or mouse marker genes may lead to inaccurate annotation. Especially, there is still a lack of special methods for economic animals (such as pigs, cows, chickens, etc.) or non-model organisms. For example, the immune cell marker genes of pigs may be different from those of humans, and existing annotation tools often lack optimization for these species. This problem limits the application of scRNA-seq in fields such as agriculture, veterinary medicine, and comparative genomics. With the expansion of single-cell data, traditional manual annotation methods have been difficult to meet the needs of high-throughput analysis. Therefore, developing more efficient and accurate cross-species single-cell annotation methods, especially special analysis processes for economic animals, has become an important direction of current single-cell research. SUMMARY

[0006] In order to overcome the problems in the prior art, the purpose of the present application is to provide an animal single-cell data cell type annotation method and system, which can accurately reflect the cell expression characteristics of the target species and improve the accuracy and reliability of cell type annotation. To achieve the above purpose, the present application provides an animal single-cell data cell type annotation method, comprising the following steps:

[0007] S1: Collecting Bulk RNA-seq data of various cell types of a specific animal, and integrating the pre-processed data into a reference dataset;

[0008] S2: Obtaining single-cell RNA-seq raw data, screening high-quality cells from the single-cell RNA-seq raw data, and then screening high-variant genes from the high-quality cells, processing the high-variant genes, extracting principal components, and reducing the dimensionality of the principal components, and performing cell clustering based on the dimensionality reduction results to obtain cell clusters;

[0009] S3: Screening high-variant genes based on the reference dataset, and taking the intersection with the single-cell high-variant gene set, then calculating the expression similarity between the cell clusters and the reference dataset based on the common high-variant genes, and determining the initial cell type of the cell clusters according to the similarity results, and iteratively optimizing the initial cell types with close scores to obtain the final cell type annotation results.

[0010] Further, the method of collecting Bulk RNA-seq data of purified cell types in each tissue and organ of a specific animal in step S1, and integrating the pre-processed data into a reference dataset:

[0011] Bulk RNA-seq data of multiple cell types are collected, all collected data are grouped according to cell type categories to obtain multiple groups of purified cell Bulk RNA-seq data of single cell types, and then the purified cell Bulk RNA-seq data are uniformly converted into TPM format, and the converted TPM data are logarithmically normalized to reduce the expression range, and integrated to form a reference data set.

[0012] Further, the method for uniformly converting the multiple groups of purified cell Bulk RNA-seq data of single cell types into TPM format is that the purified cell Bulk RNA-seq data are in the format of original gene expression Count, normalized FPKM or TPM data, and different formats are uniformly converted into TPM format.

[0013] The conversion formula for converting the original gene expression Count into TPM is as follows:

[0014]

[0015] wherein, C g is the original gene expression of gene g, l g is the exon length of gene g, i is the gene index number, is the sum of all original gene expressions in the sample divided by the exon length;

[0016] The conversion formula for converting FPKM into TPM is as follows:

[0017]

[0018] wherein, FPKM g is the FPKM standardized value of gene g, ∑ i FPKM g is the sum of FPKM values of all genes in the sample;

[0019] The formula for logarithmically normalizing the converted TPM data is as follows:

[0020]

[0021] wherein, TPM gs is the TPM expression of gene h in sample s, and sf is a scaling factor.

[0022] Further, in the step S2, the method for screening high-quality cells from the single cell RNA-seq raw data and then screening high-variant genes from the high-quality cells is as follows:

[0023] Extracting the total gene expression, the number of expressed genes and the proportion of mitochondrial genes of each cell from the single-cell RNA-seq raw data, presetting the upper threshold and the lower threshold of the total gene expression and the number of expressed genes, the lower threshold being a real number, and the upper threshold being the median plus 3 times the interquartile range, wherein the interquartile range is the difference between the upper quartile and the lower quartile; presetting the threshold of the proportion of mitochondrial genes, when the proportion of mitochondrial genes in the cell is greater than the threshold, the cell is considered unqualified; removing all low-quality cells that do not meet the threshold requirements to obtain high-quality cells;

[0024] Log-normalizing the raw gene expression of the high-quality cells to eliminate the expression deviation caused by library differences, the formula being:

[0025]

[0026] Wherein, X ij is the raw gene expression of gene j in cell i, scale_factor is the scaling factor, and ∑ k X ik is the sum of all genes k in cell i;

[0027] The high-variation genes are screened by the standard deviation method, wherein the first N genes with the most significant gene expression difference are selected.

[0028] Further, in step S2, the method for extracting principal components and reducing the dimension of the principal components after processing the high-variation genes is:

[0029] The high-variation genes are standardized, the formula being:

[0030]

[0031] Wherein, X norm，ij is the log-normalized value of gene j in cell i; μ j is the mean of gene j, and σ j is the standard deviation of gene j;

[0032] Based on the standardized high-variation gene expression data, principal component analysis is performed to reduce the dimension, generating a series of principal components and a scree plot for determining the number of principal components, and a number of principal components are selected by observing the scree plot;

[0033] Based on the selected principal components, the UMAP algorithm is used to map the cells to a two-dimensional plane, and the cells are plotted in a coordinate system with UMAP values as the horizontal and vertical coordinates.

[0034] Further, in the step S2, the method of performing cell clustering based on the dimension reduction result to obtain cell clusters is: calculating the similarity between cells based on the principal component data obtained after dimension reduction, and constructing a K-neighbor graph between cells, wherein the number K of neighbors is determined by a preset parameter;

[0035] Then, the Louvain algorithm or Leiden algorithm is used to perform iterative community division on the cells in the K-neighbor graph, the community to which the cells belong is adjusted to maximize the modularity index of the K-neighbor graph, and the division result is converged until the single-cell clustering result is obtained; and the clustering effect is verified by the two-dimensional visualization result obtained by UMAP after clustering.

[0036] Further, in the step S3, the method of filtering high-variation genes based on the reference dataset, and taking the intersection with the single-cell high-variation gene set, and then calculating the expression similarity between the cell clusters and the reference dataset based on the common high-variation genes, and determining the initial cell type of the cell cluster according to the similarity result is:

[0037] The high-variation genes are filtered based on the reference dataset by the standard deviation method or the median difference fold change method, and the intersection with the single-cell high-variation gene set is taken to obtain common high-variation genes, the expression amounts of all common high-variation genes in the cell cluster are summed to obtain the total gene expression amount of the cell cluster, and the Spearman correlation coefficient between the total gene expression amount of the cell cluster and the gene expression amount of the reference dataset is calculated based on the filtered common high-variation genes;

[0038] The reference dataset cells are grouped according to the cell type, a preset quantile of the correlation coefficient of each group is taken as the cell type score, and the cell type with the highest score is taken as the initial type of the corresponding cell cluster.

[0039] Further, the method of filtering high-variation genes by the standard deviation method is: calculating the standard deviation of the expression amount of each gene in all samples; after obtaining the standard deviation, the genes are arranged in descending order of standard deviation, and the first N genes with the largest standard deviation are taken as high-variation genes;

[0040] The median difference fold change method is to identify marker genes based on the logarithmic fold change of the median expression between cell types, specifically, each cell type is compared with other cell types, the median of the expression amount of a gene in a cell type and the remaining cell type samples is calculated, and the difference fold is calculated, and a number of genes with the largest difference fold are selected as high-variation genes of the cell type; the expression is calculated as follows:

[0041]

[0042] Wherein, K is the total number of cell types in the reference dataset, N is the number of high-variation genes to be filtered, round is the rounding function, A is the base coefficient, and B is the attenuation coefficient.

[0043] The Spearman correlation coefficient is calculated based on the ranking values of the gene expression amounts, and the formula is as follows:

[0044]

[0045] Wherein, Rank(X) is the ranking value of the total expression amount of the genes of the cell cluster, Rank(Y) is the ranking value of the expression amount of the genes of the reference data set, sigma is the standard deviation, and Cov is the covariance.

[0046] Further, in the step S3, the method for iteratively optimizing the initial cell types with close scores to obtain the final cell type annotation result is as follows:

[0047] A cell type score difference threshold is preset, when the difference between the maximum score and the second highest score in the obtained initial cell types is less than the preset cell type score difference threshold, the cell types with close scores are retained, the high-variation gene screening step in the step S2 is re-executed based on the retained cell types to obtain targeted high-variation genes, and the similarity calculation process in the step S3 is repeated based on the newly screened high-variation genes until only two cell types are left, and the difference between the maximum score and the second highest score exceeds the preset cell type score difference threshold, and the cell type with a higher score is taken as the final annotation result of the cell cluster.

[0048] Further, the method further comprises S4: based on the final cell type annotation result, difference gene analysis is performed on different cell types to identify marker genes of each cell type, and functional enrichment analysis is performed on the identified marker genes to explore the biological functions thereof.

[0049] Further, in the step S4, the identification method of the marker genes is as follows:

[0050] Based on the final cell type annotation result obtained in the step S3, the log2 fold change (log2FC) and Wilcoxon test p value of the genes between different cell types are calculated, a p value threshold and an absolute value threshold of the log2 fold change are set, and the genes that simultaneously satisfy the conditions of the p value being less than the preset threshold and the absolute value of the log2 fold change being greater than the preset threshold are determined as the marker genes of the cell type.

[0051] Further, in the step S4, the method for performing functional enrichment analysis is as follows: based on the obtained marker gene set of a specific cell type, a known biological function database is obtained for comparison, a statistical method is used to calculate the probability of the genes appearing in each type of function item, a p value reflecting the enrichment significance is obtained, and finally, the function items with a p value less than a preset threshold are screened out.

[0052] The application also provides an animal single-cell data cell type annotation system, which comprises:

[0053] a reference dataset construction module for collecting and standardizing Bulk RNA-seq data of purified cell types;

[0054] a single-cell data preprocessing module for quality control, gene expression normalization and high-variant gene screening of single-cell RNA-seq data;

[0055] a cell clustering module for dimension reduction and clustering of the preprocessed single-cell data to obtain cell clusters;

[0056] a cell annotation module for calculating the similarity of cell clusters to reference datasets based on high-variant genes to determine the initial cell type of the cell clusters, and for secondary screening and similarity calculation of cell types with close scores in the initial annotation to optimize the annotation results;

[0057] a marker gene identification module for identifying marker genes of specific cell types through differential expression analysis;

[0058] a gene function enrichment module for functional enrichment analysis of the marker genes;

[0059] an interactive visualization module for constructing a single-cell visualization system based on the R language Shiny framework to display cell type distribution, gene expression characteristics and functional enrichment results.

[0060] The present application focuses on key differential characteristics of cells and filters noise by screening high-variant genes, retains core variant information and reduces the dimension of data by principal component analysis (PCA) dimension reduction, and can cluster based on gene expression patterns to ensure that cells of the same type are clustered by constructing a cell similarity network using KNN and Louvain cell clustering. Finally, the similarity between cells is calculated based on the ranking trend based on the high-variant genes by comparing the cell clusters with the reference datasets using the Spearman correlation coefficient, and the scoring of cell types close to each other is distinguished by iterative optimization to ensure the accuracy of the final annotation. The present application not only realizes automatic analysis, but also greatly improves the accuracy and reliability of cell type annotation; BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 a cell type annotation flowchart provided for Example 1 of the present application;

[0062] Figure 2 an expression heat map of the reference dataset of chickens provided for Example 1 of the present application on known marker genes;

[0063] Figure 3 a UMAP plot of single cells provided for Example 1 of the present application, with color representing sample labels of cell sources;

[0064] Figure 4UMAP plot of single cells for the present embodiment 1, color representing cluster label;

[0065] Figure 5 Manual annotation result of single cell data according to reference dataset of chicken for the present embodiment 1;

[0066] Figure 6 Automatic annotation result of single cell data for the present embodiment 1;

[0067] Figure 7 Volcano plot of differential genes for the present embodiment 1;

[0068] Figure 8 UMAP plot of single cells for the present embodiment 1, color depth representing gene expression;

[0069] Figure 9 Gene expression violin plot for the present embodiment 1;

[0070] Figure 10 Gene expression bubble plot for the present embodiment 1;

[0071] Figure 11 Gene enrichment bubble plot for the present embodiment 1;

[0072] Figure 12 Cell information and gene information UMAP plot module in single cell visualization system for the present embodiment 2;

[0073] Figure 13 Gene co-expression UMAP plot module in single cell visualization system for the present embodiment 2;

[0074] Figure 14 Violin plot or box plot module in single cell visualization system for the present embodiment 2;

[0075] Figure 15 Cell proportion plot module in single cell visualization system for the present embodiment 2;

[0076] Figure 16 Gene expression bubble plot or heat map module in single cell visualization system for the present embodiment 2;

[0077] Figure 17 Differential gene volcano plot module in single cell visualization system for the present embodiment 2. DETAILED DESCRIPTION

[0078] Embodiment 1

[0079] The present application provides an animal single cell data cell type annotation method, which comprises the following steps:Figure 1 as shown, comprising the following steps:

[0080] S1: Collect Bulk RNA-seq data of purified cell types in each tissue and organ of a specific animal, and integrate the pre-processed data into a reference dataset.

[0081] In step S1, the method of collecting Bulk RNA-seq data of purified cell types in each tissue and organ of a specific animal, and integrating the pre-processed data into a reference dataset is as follows:

[0082] In this embodiment, the main purpose is to annotate and analyze marker genes of economic animal cells, especially to accurately annotate the cell types of a specific species. First, we need to collect the purified Bulk RNA-seq data of various cell types of the species. Bulk RNA-Seq is a method of sequencing total RNA extracted from tissues, organs, and cell populations. It obtains the average expression level of a group of cells at each gene, and is used to compare the expression differences between different individuals or different tissues of the same individual. However, for systems with strong internal cell heterogeneity, the information of gene expression will be lost. Since Bulk RNA-seq cannot distinguish the heterogeneity between cells in the sample, mixed samples usually contain multiple cell types, and the expression pattern of a certain cell type cannot be obtained. In order to ensure the accuracy of cell annotation, we group the collected Bulk RNA-seq data of multiple cell types by cell type category, and obtain multiple sets of purified cell Bulk RNA-seq data of single cell types, excluding mixed cell type data.

[0083] Convert the multiple sets of purified cell Bulk RNA-seq data of single cell types into TPM format, and perform logarithmic normalization on the converted TPM data to reduce the expression range, and integrate to form a reference dataset.

[0084] Specifically, after aligning and quantifying the data under Bulk RNA-seq, we will get the expression of each gene (Count) for each sample. Generally, standardization will be performed before subsequent analysis. However, the standardization methods used in RNA-seq collected from public databases are not uniform. The most common ones are FPKM (Fragments Per Kilobase Million) or TPM (Transcripts Per Million). Therefore, we need to convert the purified cell Bulk RNA-seq data format into raw gene expression Count, standardized FPKM or TPM data, and all convert them into TPM format.

[0085] The conversion formula for converting raw gene expression Count to TPM is as follows:

[0086]

[0087] wherein, C g is the original gene expression of gene g, l g is the exon length of gene g, ∑ i C i is the total number of sequencing fragments of the sample, i is the number of gene index, is the sum of the original gene expression of all genes in the sample divided by the exon length.

[0088] The conversion formula for converting the original data FPKM to TPM is as follows:

[0089]

[0090] wherein, FPKM g is the FPKM standardized value of gene g, ∑ i FPKM g is the sum of the FPKM values of all genes in the sample.

[0091] The original data Count can also be converted to FPKM and then to TPM. The conversion formula for converting Count to FPKM is as follows:

[0092]

[0093] wherein, C g is the original gene expression Count of gene g, l g is the exon length (unit: kilobase pair kbp) of gene g, ∑ i C i is the total number of sequencing fragments of the sample (the sum of Count of all genes).

[0094] The formula for logarithmic normalization of the converted TPM data is:

[0095]

[0096] wherein, TPM gs represents the TPM expression of gene h in sample s, and sf is the scaling factor.

[0097] In this embodiment, taking the construction of the reference dataset of chicken as an example, after the reference dataset is constructed according to the above steps, a heat map is used to show the expression of the reference dataset on the known marker genes, as shown in Figure 2 The results show that each cell type of the reference dataset is specifically highly expressed on its marker genes.

[0098] S2: obtaining single-cell RNA-seq raw data, screening high-quality cells from the single-cell RNA-seq raw data, screening high-variation genes from the high-quality cells, extracting principal components after processing the high-variation genes, and performing dimension reduction on the principal components to obtain cell clusters based on the dimension reduction results;

[0099] The method for screening high-quality cells from the single-cell RNA-seq data and screening high-variation genes from the high-quality cells comprises the following steps:

[0100] The total gene expression amount (nCount), the number of expressed genes (nFeature), and the mitochondrial gene proportion (percent.mito) of each cell are extracted from the single-cell RNA-seq raw data, preset upper and lower threshold values of the total gene expression amount and the number of expressed genes, and the lower threshold value is 200 in this embodiment, and low-quality cell filtering is performed; since data from different experiments or tissues can cause significant differences in gene expression, a fixed value cannot be simply used to filter cells, therefore, the upper threshold value is flexibly set, specifically, the upper threshold value is the median value plus 3 times the interquartile range, wherein the interquartile range is the difference between the upper quartile and the lower quartile, and the median value plus 3 times the interquartile range can be considered as a larger outlier in the data. A mitochondrial gene proportion threshold value is preset, and the mitochondrial gene proportion threshold value is set to 30% in this embodiment, when the mitochondrial gene proportion in a cell is greater than 30%, the cell is considered unqualified, and low-quality cells that do not meet the threshold value requirement are removed to obtain high-quality cells.

[0101] The high-quality cells are subjected to logarithmic normalization processing to eliminate expression deviation caused by library difference, and the formula is represented as:

[0102]

[0103] wherein, X ij is the original count of the original gene j in the cell i, scale_factor is a scaling factor, which is 10000 by default in this embodiment, and ∑ k X ik is the sum of all genes k in the cell i.

[0104] In order to calculate the similarity between the single-cell gene expression profile of an economic animal and the expression profile of a pure cell type sample, generally, all genes are not used for similarity calculation, but high-variation genes are screened out first. These genes can best reflect the differences between different cell types, which can improve the accuracy of cell annotation on the one hand, and reduce the amount of calculation on the other hand. There are various ways to screen high-variation genes. High-variation genes can be screened from the gene expression data of a reference dataset, which has a large difference in expression between samples, or high-variation genes can be selected from single-cell expression data.

[0105] In this embodiment, the standard deviation method is used to screen high-variation genes (HVG) for high-quality cells. The high-variation genes contain the main difference information between cells.

[0106] In this embodiment, the preferred base coefficient is 500, and the decay coefficient is 2 / 3. The base resolution is ensured by the minimum threshold value of 500, the number of cell types is dynamically adapted by the decay coefficient of 2 / 3, the result is ensured to be convenient for subsequent operation by rounding, and finally the goal of accurately distinguishing all cell types by using a reasonable number of marker genes is achieved.

[0107] In the step S2, the method for extracting principal components and reducing the dimension of the principal components after processing the high-variation genes is:

[0108] The high-variation genes are standardized, and the formula is:

[0109]

[0110] Wherein, X norm，ij is the logarithmic normalized value of gene j in cell i; μ j is the mean of gene j, and σ j is the standard deviation of gene j.

[0111] Based on the standardized high-variation gene expression data, principal component analysis is performed for dimension reduction, a series of principal components and a gravel map for judging the number of principal components are generated, and a number of principal components are selected by observation;

[0112] Based on the selected principal components, the UMAP algorithm is used to map the cells to a two-dimensional plane, and the cells are plotted in a coordinate system with UMAP values as horizontal and vertical coordinates.

[0113] In step S2, the method for obtaining cell clusters by clustering cells based on the dimension reduction result is:

[0114] The similarity between cells is calculated based on the principal component data obtained after dimension reduction, and a K-neighbor graph between cells is constructed, wherein the number K of neighbors is determined by a preset parameter, and in this embodiment, the default is 20. Specifically, the Euclidean distance of cells in the low-dimensional space is calculated, 20 nearest neighbors of each cell are connected to form an initial neighborhood relationship; the connection between cells is weighted using Jaccard similarity and other metrics, and edges that significantly share neighbors are retained to enhance robustness.

[0115] The Louvain algorithm or Leiden algorithm is further used to perform iterative community division on the cells in the K-neighbor graph, and the modularity index of the K-neighbor graph is maximized by adjusting the community to which the cells belong until the division result converges, to obtain a single-cell clustering result. The convergence condition of the iterative community division is that the change amount of the modularity in two consecutive iterations is less than a preset threshold, for example, the threshold is set to 1x10 -5 -4, or the number of iterations reaches a preset upper limit of 100 times. Finally, a cell cluster with biological significance is output.

[0116] After clustering, the two-dimensional visualization result obtained by UMAP is used to verify the clustering effect. The UMAP visualization verifies the spatial separation, and the whole process integrates local distance measurement and global graph structure analysis, balancing the calculation efficiency and clustering accuracy. As shown in FIG. 6, the horizontal and vertical coordinates are UMAP values, each point in the figure represents a cell, the distance between cells reflects the similarity of gene expression, and the color of each point represents the sample from which the cell is derived. As shown in FIG. 7, the color of the point represents the cell cluster label after clustering. It can be seen that the cells in the same cell cluster are closer to each other, and the cells in different cell clusters are farther away from each other. Figure 3 Figure 4

[0117] S3: filtering high-variant genes based on the reference dataset, taking intersection with the single-cell high-variant gene set, then calculating the expression similarity between the cell cluster and the reference dataset based on the common high-variant genes, determining the initial cell type of the cell cluster according to the similarity result, and iteratively optimizing the initial cell types with close scores to obtain the final cell type annotation result.

[0118] The principle of cell type annotation is that cells of the same type have higher similarity in gene expression, and using a reference dataset for annotation is to find which sample expression profile of a cell type is most similar to the single-cell expression profile, and the cell is annotated as the cell type.

[0119] In step S3, the method of filtering high-variant genes based on the reference dataset, taking intersection with the single-cell high-variant gene set, then calculating the expression similarity between the cell cluster and the reference dataset based on the common high-variant genes, and determining the initial cell type of the cell cluster according to the similarity result is:​​

[0120] Filter high variability genes based on reference dataset by standard deviation method or median fold change method, and get common high variability genes by intersecting with single cell high variability gene set.

[0121] The method using standard deviation method is: calculate the standard deviation of the expression of each gene in all samples; after obtaining the standard deviation, arrange the genes in descending order of standard deviation, and take the first N genes with the largest standard deviation as high variability genes. This can better distinguish different cell types and is conducive to accurate annotation of cell types. After obtaining the standard deviation, arrange the genes in descending order of standard deviation, and the default setting of N is 2000.

[0122] The median fold change method is based on the log-fold change of the median expression between cell types to identify marker genes for screening. Specifically, each cell type is compared with other cell types, the median of the expression of the gene in a certain cell type and the remaining cell type samples is calculated, and the fold change is calculated. Select the largest number of genes as the high variability genes of the cell type; the number of high variability genes can be automatically set according to the number of cell types, and the gene number calculation expression is as follows:

[0123]

[0124] Wherein, K is the total number of cell types in the reference dataset, N is the number of high variability genes to be screened, round is the rounding function, A is the base coefficient, and B is the attenuation coefficient.

[0125] Sum the expression of all common high variability genes in the cell cluster to get the total gene expression of the cell cluster, and calculate the Spearman correlation coefficient between the total gene expression of the cell cluster and the gene expression of the reference dataset based on the screened common high variability genes. The Spearman correlation coefficient is a non-parametric rank correlation coefficient, which is used to measure the monotonic relationship between two variables. It is calculated based on the ordering values of the variables rather than the original data values, so it is not sensitive to outliers and is suitable for nonlinear but monotonic relationship analysis. The present application calculates the similarity between single cell and sample by means of Spearman correlation coefficient. When the expression profiles of the cell and the sample have consistent monotonicity, it means that the cell and the sample have higher consistency in gene expression, that is, they simultaneously highly express and simultaneously lowly express in some genes. Therefore, the Spearman correlation coefficient can effectively measure the similarity between the expression profiles of the cell and the sample.

[0126] When calculating the Spearman correlation coefficient, the high variability gene expression of the cell and the sample needs to be sorted from small to large respectively, and the ordering value Rank(X) is used instead of the original expression value X. The formula of the Spearman correlation coefficient is:

[0127]

[0128] wherein Rank(X) is the rank value of the total expression of the genes in the cell cluster, Rank(Y) is the rank value of the expression of the genes in the reference dataset, σ is the standard deviation, and Cov is the covariance.

[0129] According to the cell type grouping of the reference dataset, a preset quantile of the correlation coefficient of each group is taken as the cell type score, and the cell type with the highest score is taken as the initial type of the corresponding cell cluster. In this embodiment, the preset quantile of the correlation coefficient of each group is set to 80%.

[0130] The above method can obtain the score of each cell cluster of single-cell data on each cell type, however, in practice, the scores of several cell types with the highest scores may have a small difference, which can cause the problem of inaccurate cell type annotation, and therefore further iteration and optimization of the scores of several cell types with close scores is needed.

[0131] In the step S3, the method for iteration and optimization of the initial cell types with close scores to obtain the final cell type annotation result is as follows: a preset cell type score difference threshold is set, and in this embodiment, the cell type score difference threshold is set to 0.05; when the difference between the maximum score and the second highest score in the obtained initial cell types is less than 0.05, the cell types with close scores are retained, the high-variation gene screening step in the step S2 is re-executed based on the retained cell types to obtain targeted high-variation genes, and the similarity calculation process in the step S3 is repeated based on the newly screened high-variation genes until only two cell types are left, and the difference between the maximum score and the second highest score exceeds the preset cell type score difference threshold 0.05, and the cell type with the higher score is taken as the final annotation result of the cell cluster.

[0132] Taking the single-cell dataset CRA002353 of chicken as an example, first, the dataset is manually annotated using known marker genes, and the annotation result is as shown in Figure 5 , and then the dataset is automatically annotated using the method based on the chicken reference dataset constructed in the step S1, and the annotation result is as shown in Figure 6 From the figure, it can be seen that the automatic annotation result is consistent with the manual annotation result on the main cell types.

[0133] S4: Based on the final cell type annotation result, differential gene analysis is performed on different cell types to identify the marker genes of each cell type, and functional enrichment analysis is performed on the identified marker genes to explore their biological functions.

[0134] In step S4, the identification of the marker gene includes: based on the final cell type annotation results obtained in step S3, calculating the logarithmic fold change (log2FC) and the Wilcoxon p-value between different cell types. The method for calculating the logarithmic fold change (log2FC) is as follows: for each gene, calculate the average expression level of the gene in any two cell clusters, then calculate the two averages to obtain the fold change (FC), and then take the logarithm of the fold change to obtain the logarithmic fold change (log2FC). When this value is positive, it indicates that the gene expression level in the first cluster is upregulated compared to the second cluster; conversely, when this value is negative, it indicates that the gene expression level in the first cluster is downregulated compared to the second cluster.

[0135] By setting thresholds for the p-value and the absolute value of the logarithmic fold change, genes that simultaneously satisfy the condition of a p-value less than the preset threshold and an absolute value of the logarithmic fold change greater than the preset threshold are identified as marker genes for that cell type. In this embodiment, a p-value less than 0.05 is considered significant, and differentially expressed genes are selected based on the condition |log2FC|>1.

[0136] Taking the dataset CRA002353 as an example, calculate the differentially expressed genes in myoblasts compared to all other cells, such as... Figure 7 As shown, the horizontal axis is log2FC and the vertical axis is -log 10 p value The dashed line represents the threshold for screening differentially expressed genes. The red part on the right represents the upregulated genes in myoblasts, and the blue part on the left represents the downregulated genes in myoblasts.

[0137] Furthermore, the expression characteristics of marker genes are visualized. Taking the RBM24 gene in myoblasts as an example, the expression level of this gene is mapped to the colors of a UMAP plot, such as... Figure 8 As shown, RBM24 is the darkest color in the myoblast cluster, indicating that this gene is specifically highly expressed only in myoblasts. Figure 9 The diagram shown is a violin plot of gene expression. The distribution of RBM24 gene expression in myoblasts is wider at the top and narrower at the bottom, indicating that most myoblasts highly express the RBM24 gene. Figure 10 This is a bubble graph of gene expression. The color depth of the bubbles represents the average gene expression level in the cell cluster, and the bubble size represents the proportion of cells in the cluster that express that gene. The graph shows that MYOD1, RBM24, and TMSB15B genes are specifically highly expressed in myoblasts. The above statistical graphs confirm that RBM24 is a marker gene for myoblasts.

[0138] After identifying the single-cell marker genes, further functional exploration is essential. Gene enrichment analysis is a statistical method for identifying whether a set of genes is significantly over-represented in a particular biological function, pathway or regulatory network. Through gene enrichment analysis, the biological processes, molecular functions and signal pathways involved in the marker genes can be revealed, thereby deepening the understanding of their potential regulatory mechanisms in cells.

[0139] The method for performing functional enrichment analysis is as follows: based on the obtained marker gene set of a specific cell type, a known biological function database is obtained for comparison. The biological function database, such as the GO or KEGG database, calculates the probability of the occurrence of these genes in each functional item using statistical methods, obtains the p-value reflecting the enrichment significance, and finally screens out the functional items with a p-value less than a preset threshold.

[0140] Specifically, gene enrichment analysis of the marker genes of a certain cell type can explore the active biological pathways in the cell. Again taking myoblasts as an example, the marker gene set is screened out, and the R language clusterProfiler package enrichGO function is used to perform GO enrichment analysis on the set, to obtain the number of genes, the proportion of genes, the functional items and the significance p-value. The enrichment result is shown in Figure 11 The vertical axis shows the most significantly enriched pathway items, the horizontal axis represents the proportion of enriched genes, the bubble size represents the number of enriched genes, the bubble color represents the p-value, and the enriched pathways are sorted according to the proportion of enriched genes and the significance p-value. It can be seen that the marker genes of myoblasts are mainly enriched in the pathways related to skeletal muscle development and differentiation.

[0141] Example 2

[0142] The present application also provides an animal single-cell data cell type annotation system, which is applied to the animal single-cell data cell type annotation method described in the above example 1, and comprises:

[0143] A reference dataset construction module is configured to collect Bulk RNA-seq data of purified cell types and perform standardization processing;

[0144] A single-cell data preprocessing module is configured to perform quality control, gene expression normalization and high-variation gene screening on single-cell RNA-seq data;

[0145] A cell clustering module is configured to perform dimension reduction and clustering on the preprocessed single-cell data to obtain cell clusters;

[0146] A cell annotation module is configured to calculate the similarity of cell clusters and reference datasets based on high-variation genes, determine the initial cell type of the cell clusters, perform secondary screening and similarity calculation on the cell types with close scores in the initial annotation, and optimize the annotation result.

[0147] Marker gene identification module, for identifying marker genes of specific cell types through differential expression analysis;

[0148] Gene function enrichment module, for functional enrichment analysis of marker genes;

[0149] Interactive visualization module, based on the R language Shiny framework to build a single-cell visualization system, for displaying cell type distribution, gene expression characteristics and functional enrichment results.

[0150] The interactive visualization module is based on the above analysis results to build a single-cell visualization system, which is based on the R language Shiny framework, and is convenient for researchers to view the analysis results of economic animal single-cell data. The system is characterized by interactive viewing of single-cell data visualization results, easy viewing of interested genes or cell types, and control of statistical chart elements such as size and color. The interactive visualization module is divided into 7 sub-modules, which are cell information UMAP chart, gene information UMAP chart, gene co-expression module, violin chart or box plot module, cell proportion module, bubble chart or heat map module, and differential gene volcano plot module. The cell information UMAP chart is shown as Figure 12 The left side, mainly showing the classification label information of the cells or the continuous indicators about the cells, the cell label distribution or the distribution of different cells on the continuous indicators can be easily viewed, and the cell-related variable name can be selected in the drop-down box on the left side of the figure, and the multiple selection box can select to display part of the cell clusters. The gene information UMAP chart is shown as Figure 12 The right side, which can intuitively view the expression amount of each cluster of cells on the interested gene in the UMAP coordinate system, and the interested gene name can be input or selected in the drop-down box above the figure. The gene co-expression module is shown as Figure 13 The distribution of the expression amounts of two genes can be displayed in one UMAP chart, and two different gene names can be input in the two drop-down boxes on the left side. The violin chart or box plot module is shown as Figure 14 The module displays the distribution of gene expression or cell continuous indicators in different cell clusters, the first drop-down box on the left side can select different cell classification variables to change the classification name of the horizontal coordinate, the second drop-down box on the left side can select the continuous variable about the cells or the interested gene name, and the two options on the left side can control the display of the violin chart and the box plot. The cell proportion module is shown as Figure 15 The proportion of different cell clusters is shown, and different cell classification labels can be selected in the two drop-down boxes on the left side, and the two options on the left side can control the display of the cell proportion or the absolute number of cells. The bubble chart or heat map module is shown as Figure 16As shown, the module can simultaneously show the expression of multiple genes in different cell clusters. The volcano plot module is shown in FIG. 6B. Figure 17 As shown, the module shows the differentially expressed genes for each cell cluster, and the different cell classification variables, cell label names, log fold change, and significance p-value thresholds can be selected in the left drop-down box.

Claims

1. A method for annotating cell types in animal single-cell data, characterized in that, Includes the following steps: S1: Collect Bulk RNA-seq data of multiple cell types in specific animals, and integrate them into a reference dataset after preprocessing; S2: Obtain raw single-cell RNA-seq data, screen the raw single-cell RNA-seq data to obtain high-quality cells, then screen out highly variable genes from the high-quality cells, process the highly variable genes, extract principal components and reduce the dimensionality of the principal components, and perform cell clustering based on the dimensionality reduction results to obtain cell clusters. S3: Based on the reference dataset, highly variable genes are screened and their intersection with the single-cell highly variable gene set is taken. Then, the expression similarity between the cell cluster and the reference dataset is calculated based on the common highly variable genes. The initial cell type of the cell cluster is determined according to the similarity results. The initial cell types with similar scores are iteratively optimized to obtain the final cell type annotation results.

2. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, The method described in step S1 for collecting Bulk RNA-seq data of purified cell types from various tissues and organs of specific animals, and integrating them into a reference dataset after preprocessing: Bulk RNA-seq data of multiple cell types were collected, and all collected data were grouped according to cell type category to obtain multiple groups of purified cell Bulk RNA-seq data of a single cell type; Then, the purified cell Bulk RNA-seq data were uniformly converted into TPM format, and the converted TPM data were log-normalized to narrow the expression level range, and integrated to form a reference dataset.

3. The method for annotating cell types in animal single-cell data according to claim 2, characterized in that, The method for uniformly converting multiple sets of purified cell Bulk RNA-seq data of a single cell type into TPM format is as follows: the purified cell Bulk RNA-seq data format is the original gene expression level Count, the standardized FPKM or TPM data, and the different formats are uniformly converted into TPM format; The conversion formula for converting the original gene expression level Count to TPM is as follows: Among them, C g The original gene expression level of gene g, l g Let be the exon length of gene g, and i be the gene index number. The sum of the expression levels of all original genes in the sample, excluding the length of the original genes; The conversion formula for FPKM to TPM is as follows: Among them, FPKM g Let ∑ be the FPKM normalized value of gene g. i FPKM g This is the sum of the FPKM values ​​of all genes in the sample; The formula for logarithmic normalization of the converted TPM data is expressed as follows: Among them, TPM gs To represent the TPM expression level of gene h in sample s, sf is a scaling factor.

4. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, In step S2, the method of screening the raw single-cell RNA-seq data to obtain high-quality cells, and then screening for highly variable genes from the high-quality cells, is as follows: From the raw single-cell RNA-seq data, the total gene expression level, the number of expressed genes, and the proportion of mitochondrial genes in each cell were extracted. Upper and lower thresholds were preset for the total gene expression level and the number of expressed genes. The lower threshold was a real number, and the upper threshold was the median plus three times the interquartile range, where the interquartile range was the difference between the upper and lower quartiles. A threshold was also preset for the mitochondrial gene proportion; cells with a mitochondrial gene proportion greater than this threshold were considered unqualified. All low-quality cells that did not meet the threshold requirements were removed, resulting in high-quality cells. Log-normalization was performed on the original gene expression levels of high-quality cells to eliminate expression level bias caused by library differences. The formula is as follows: Among them, X ij ∑ represents the original gene expression level of gene j in cell i, where scale_factor is the scaling factor. k X ik Summing the sum of all genes k within cell i; High-variance genes were screened using the standard deviation method for high-quality cells, where the standard deviation method selected the top N genes with the most significant differences in gene expression levels.

5. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, In step S2, the method for processing highly variable genes, extracting principal components, and then reducing the dimensionality of the principal components is as follows: The highly variable genes are standardized using the following formula: Among them, X norm,ij It is the log-normalized value of gene j in cell i; μ j Let σ be the mean of gene j. j The standard deviation of gene j; Principal component analysis was performed to reduce the dimensionality of the standardized high-variability gene expression data, generating a series of principal components and a scree plot to determine the number of principal components. Several principal components were then selected by observing the scree plot. Based on the selected principal components, the UMAP algorithm is used to map the cells to a two-dimensional plane, and the cells are drawn in a coordinate system with UMAP values ​​as the horizontal and vertical coordinates.

6. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, In step S2, the method for obtaining cell clusters by cell clustering based on the dimensionality reduction results is as follows: The similarity between cells is calculated based on the principal component data obtained after dimensionality reduction, and a K-nearest neighbor graph between cells is constructed, where the number of nearest neighbors K is determined by preset parameters; Then, the Louvain algorithm or Leiden algorithm is used to iteratively divide the cells in the K-nearest neighbor graph into communities. By adjusting the communities to which the cells belong, the modularity index of the K-nearest neighbor graph is maximized until the division results converge, and the single-cell clustering results are obtained. After clustering is completed, the clustering effect is verified by two-dimensional visualization results obtained through UMAP.

7. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, In step S3, the method of screening highly variable genes based on a reference dataset, taking the intersection with the single-cell highly variable gene set, calculating the expression similarity between the cell cluster and the reference dataset based on the common highly variable genes, and determining the initial cell type of the cell cluster based on the similarity results is as follows: Highly variable genes were screened using the standard deviation method or the median difference fold method based on the reference dataset. Common highly variable genes were obtained by taking the intersection with the single-cell highly variable gene set. The expression levels of all common highly variable genes in the cell cluster were summed to obtain the total gene expression level of the cell cluster. Based on the screened common highly variable genes, the Spearman correlation coefficient between the total gene expression level of the cell cluster and the gene expression level of the reference dataset was calculated. The cells are grouped according to the cell types in the reference dataset. The preset quantile of the correlation coefficient of each group is taken as the cell type score, and the cell type with the highest score is taken as the initial type of the corresponding cell cluster.

8. The method for annotating cell types in animal single-cell data according to claim 7, characterized in that, The method for screening highly variable genes using the standard deviation method is as follows: calculate the standard deviation of the expression level of each gene in the sample; after obtaining the standard deviation, sort the genes in descending order of standard deviation, and take the top N genes with the largest standard deviation as highly variable genes. The median fold change method is used to identify marker genes based on the logarithmic fold change in median expression among cell types. Specifically, each cell type is compared with other cell types, the median expression level of the gene in a particular cell type and the remaining cell types is calculated, and then the fold change is calculated. The genes with the largest fold changes are selected as highly variable genes for that cell type. The calculation expression is as follows: Where K is the total number of cell types in the reference dataset, N is the number of highly variable genes to be screened, round is the rounding function, A is the base coefficient, and B is the decay coefficient; The Spearman correlation coefficient is calculated based on the gene expression level ranking, and the formula is as follows: Where Rank(X) is the ranking value of the total gene expression level of the cell cluster, Rank(Y) is the ranking value of the gene expression level of the reference dataset, σ is the standard deviation, and Cov is the covariance.

9. The method for annotating cell types in animal single-cell data according to claim 8, characterized in that, In step S3, the method for iteratively optimizing the initial cell types with similar scores to obtain the final cell type annotation results is as follows: A preset cell type difference threshold is set. When the difference between the highest and second-highest scores in the initial cell types is less than the preset cell type difference threshold, these cell types with similar scores are retained. Based on the retained cell types, the high-variable gene screening step in step S2 is re-executed. Then, based on the newly screened high-variable genes, the similarity calculation process in step S3 is repeated until only two cell types remain, and the difference between the highest and second-highest scores exceeds the preset cell type difference threshold. The cell type with the higher score is taken as the final annotation result of the cell cluster.

10. The method for annotating cell types in animal single-cell data according to claim 1, characterized in that, It also includes S4: Based on the final cell type annotation results, differential gene analysis is performed on different cell types to identify marker genes for each cell type, and functional enrichment analysis is performed on the identified marker genes.

11. The method for annotating cell types in animal single-cell data according to claim 10, characterized in that, In step S4, the method for identifying the marker gene is as follows: Based on the final cell type annotation results obtained in step S3, the logarithmic fold change (log2FC) and Wilcoxon test p-value of genes between different cell types are calculated. A p-value threshold and an absolute value threshold for the logarithmic fold change are set. Genes that simultaneously satisfy the condition that the p-value is less than the preset threshold and the absolute value of the logarithmic fold change is greater than the preset threshold are identified as marker genes of that cell type.

12. The method for annotating cell types in animal single-cell data according to claim 10, characterized in that, In step S4, the method for performing functional enrichment analysis is as follows: based on the obtained set of marker genes for a specific cell type, a known biological function database is obtained for comparison, statistical methods are used to calculate the probability of these genes appearing in various functional items, a p-value reflecting the significance of enrichment is obtained, and finally functional items with p-values ​​less than a preset threshold are selected.

13. A cell type annotation system for animal single-cell data, applied to the cell type annotation method for animal single-cell data as described in any one of claims 1-12, characterized in that, include: A reference dataset building module is used to collect and standardize Bulk RNA-seq data from purified cell types. The single-cell data preprocessing module is used for quality control, gene expression level normalization, and screening of highly variable genes in single-cell RNA-seq data. The cell clustering module is used to perform dimensionality reduction and clustering on preprocessed single-cell data to obtain cell clusters. The cell annotation module is used to calculate the similarity between cell clusters and reference datasets based on highly variable genes, and to determine the initial cell type of the cell clusters. A secondary screening and similarity calculation were performed on cell types with similar scores in the initial annotation to optimize the annotation results; The marker gene identification module is used to identify marker genes for specific cell types through differential expression analysis. The gene function enrichment module is used to perform functional enrichment analysis on marker genes. The interactive visualization module is a single-cell visualization system built on the Shiny framework of the R language, used to display cell type distribution, gene expression characteristics, and functional enrichment results.

Images