Cell analysis method, device and equipment for bulk data

Abstract

The embodiment of the invention relates to the technical field of bioinformatics, and provides a bulk data cell analysis method, device and equipment, and the method comprises the following steps: constructing an initial reference matrix according to a single cell data set and a cell type annotation template, each element in the initial reference matrix represents the gene expression quantity of each cell state under each characteristic gene; performing deconvolution on the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result; updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix; and according to the first reference matrix, performing deconvolution on the bulk data to be analyzed to obtain a second proportion and a second gene expression quantity of each cell type in the bulk data to be analyzed. According to the embodiment of the invention, the accuracy of cell analysis in bulk data can be improved.

Description

TECHNICAL FIELD

[0001] Embodiments of the present specification relate to the technical field of bioinformatics, and particularly to a bulk data cell analysis method, device and equipment. BACKGROUND

[0002] Tumor microenvironment (TME) has an important influence on the development of cancer. According to the characteristics and changes of TME, a personalized treatment plan can be developed to improve treatment effect and patient survival rate. TME includes various immune cell types, cancer-associated fibroblasts (CAFs), endothelial cells, pericytes and other tissue-resident cell types. Current analysis of TME usually requires invasive surgery or biopsy, which has certain risks and often leads to information lag. In recent years, deconvolution algorithms based on bulk data can improve the efficiency of microenvironment analysis to some extent, but there are still certain limitations.

[0003] Current deconvolution methods mainly include reference-based methods and unsupervised methods. Reference-based methods deconvolve bulk data by using a fixed reference matrix, for example, CIBERSORT uses 22 cell types as reference. Reference-based methods based on reference matrix make the results more interpretable and reliable, but due to incomplete cell coverage, batch effect and tumor heterogeneity, etc., the effect of reference-based deconvolution method is poor. Unsupervised methods directly identify cell types from data through clustering or mathematical models without relying on any reference to deconvolve bulk data, for example, using mathematical model mixed convex analysis to directly identify subclone marker genes from bulk gene expression. Although this unsupervised method does not require a reference matrix, it is difficult to distinguish between similar cell types and there is a problem of incomplete cell types. SUMMARY

[0004] In view of the above problems of the prior art, the purpose of the embodiments of the present specification is to provide a bulk data cell analysis method, device and equipment to improve the accuracy of cell analysis in bulk data.

[0005] To solve the above technical problems, the specific technical solutions of the embodiments of the present specification are as follows:

[0006] In one aspect, the embodiments of the present specification provide a bulk data cell analysis method, which comprises:

[0007] Obtaining a single-cell data set and a cell type annotation template, the cell type annotation template comprising a pre-defined cell type hierarchy and a feature gene set;

[0008] constructing an initial reference matrix according to the single-cell data set and the cell type annotation template, each element in the initial reference matrix representing a gene expression amount of each cell state under each characteristic gene;

[0009] deconvolving the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result, the first deconvolution result including a first proportion and a first gene expression amount of each cell type;

[0010] updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix;

[0011] deconvolving the bulk data to be analyzed according to the first reference matrix to obtain a second deconvolution result, the second deconvolution result including a second proportion and a second gene expression amount of each cell type.

[0012] Further, the constructing an initial reference matrix according to the single-cell data set and the cell type annotation template comprises:

[0013] annotating each single cell in the single-cell data set by using the cell type annotation template to obtain a cell state label of each single cell;

[0014] summing expression data of all single cells of the same cell state according to the cell state label of each single cell to obtain expression counts of each cell state under each characteristic gene;

[0015] normalizing the expression counts of each cell state under each characteristic gene to obtain gene expression amounts of each cell state under each characteristic gene;

[0016] constructing the initial reference matrix according to the gene expression amounts of each cell state under each characteristic gene.

[0017] Further, the deconvolving the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result comprises:

[0018] S1: taking the initial reference matrix and the bulk data to be analyzed as a prior distribution and sample information respectively, and determining a joint posterior distribution of cells in the bulk data to be analyzed by using a Gibbs sampling method;

[0019] S2: marginalizing the joint posterior distribution to obtain a gene expression marginal distribution matrix and a proportion marginal distribution of cell states in the bulk data to be analyzed;

[0020] S3: calculating a gene expression mean matrix and a proportion mean of cell states according to the gene expression marginal distribution matrix and the proportion marginal distribution of the cell states.

[0021] S4: summing the gene expression mean matrix and the proportion mean of all cell states belonging to the same cell type to obtain a gene expression matrix and a proportion of the cell type;

[0022] S5: updating the initial reference matrix according to the gene expression matrix of the cell type to obtain a fourth reference matrix;

[0023] S6: taking the fourth reference matrix as a new prior distribution, iteratively performing steps S1-S5 until a preset iteration termination condition is reached to obtain a first deconvolution result.

[0024] Further, the cell types include malignant cells and non-malignant cells, and the updating of the initial reference matrix according to the gene expression matrix of the cell type to obtain a fourth reference matrix includes:

[0025] According to the gene expression matrix of the cell type, a second reference matrix of the malignant cells and an initial reference matrix of the non-malignant cells are calculated;

[0026] According to the initial reference matrix of the non-malignant cells and a preset prior probability, a third reference matrix of the non-malignant cells is calculated;

[0027] The second reference matrix and the third reference matrix are merged to obtain the fourth reference matrix.

[0028] Further, the calculation of the second reference matrix of the malignant cells according to the gene expression matrix of the cell type includes:

[0029] According to the gene expression matrix of the cell type, the expression value of each feature gene of the malignant cells is calculated by using the following formula:

[0030]

[0031] wherein, represents the second reference matrix of the malignant cells in the bulk data to be analyzed, t1 represents the malignant cells, represents the gene expression matrix of the malignant cells in the bulk data to be analyzed;

[0032] The second reference matrix is constructed according to the expression value of each feature gene of the malignant cells.

[0033] Further, the calculation of the initial reference matrix of the non-malignant cells according to the gene expression matrix of the cell type includes:

[0034] According to the gene expression matrix of the cell type, the matrix value in the initial reference matrix of the non-malignant cells is calculated by using the following formula:

[0035]

[0036] wherein, represents a matrix value in the initial reference matrix of non-malignant cells, t2 represents non-malignant cells, W cg represents the original count of the c th cell in the g th gene in the non-malignant cells;

[0037] The initial reference matrix is constructed according to the matrix values in the initial reference matrix.

[0038] Further, the third reference matrix of non-malignant cells is calculated according to the initial reference matrix of non-malignant cells and the preset prior probability, and the calculation includes:

[0039] The matrix value in the reference matrix of each characteristic gene of non-malignant cells is calculated by using the following formula:

[0040]

[0041] wherein, represents a matrix value in the reference matrix of each characteristic gene of non-malignant cells, represents a matrix value in the initial reference matrix of non-malignant cells, represents a prior probability;

[0042] The reference matrix of each characteristic gene of non-malignant cells is constructed according to the matrix values in the reference matrix of each characteristic gene of non-malignant cells.

[0043] The reference matrix of each characteristic gene of non-malignant cells is merged to obtain the third reference matrix of non-malignant cells.

[0044] Further, the initial reference matrix is updated according to the first deconvolution result to obtain the first reference matrix, and the updating includes:

[0045] The cell types with a first proportion less than a first preset threshold in the first deconvolution result are screened out to obtain a first cell type set;

[0046] The cell types with a correlation less than a second preset threshold with a cancer type in the first cell type set are screened out to obtain a second cell type set;

[0047] The first reference matrix is obtained according to the first gene expression amount of each cell type in the second cell type set.

[0048] On the other hand, the embodiment of the present specification also provides a cell analysis device of bulk data, and the device includes:

[0049] The acquisition module is configured to acquire a single-cell data set and a cell type annotation template, the cell type annotation template including a predefined cell type hierarchy and a characteristic gene set.

[0050] The construction module is configured to construct an initial reference matrix according to the single-cell data set and the cell type annotation template, each element in the initial reference matrix representing a gene expression amount of each cell state under each characteristic gene.

[0051] The first deconvolution module is configured to perform deconvolution on the bulk data to be analyzed according to the initial reference matrix, to obtain a first deconvolution result, the first deconvolution result including a first proportion and a first gene expression amount of each cell type.

[0052] The update module is configured to update the initial reference matrix according to the first deconvolution result, to obtain a first reference matrix.

[0053] The second deconvolution module is configured to perform deconvolution on the bulk data to be analyzed according to the first reference matrix, to obtain a second deconvolution result, the second deconvolution result including a second proportion and a second gene expression amount of each cell type.

[0054] In yet another aspect, the embodiments of the present specification also provide a computer device, including a memory, a processor, and a computer program stored in the memory, when the computer program is run by the processor, instructions of any one of the above methods are executed.

[0055] In yet another aspect, the embodiments of the present specification also provide a computer-readable storage medium, having a computer program stored thereon, when the computer program is run by a processor of a computer device, instructions of any one of the above methods are executed.

[0056] In yet another aspect, the embodiments of the present specification also provide a computer program product, when the computer program product is run by a processor of a computer device, instructions of any one of the above methods are executed.

[0057] By adopting the technical scheme, the bulk data cell analysis method provided in the embodiments of the present specification first obtains an initial reference matrix containing multiple cell types through a single cell data set and a cell type annotation template, the initial reference matrix is composed of a tumor microenvironment (TME) cell hierarchy of a general cancer, each parent node in the cell hierarchy represents a group of related cells, and each terminal node represents a target single cell category. Then, the bulk data to be analyzed is deconvoluted based on a deconvolution algorithm and the initial reference matrix to obtain a first deconvolution result, and then an updated deconvolution reference matrix is obtained according to the first deconvolution result. The updated deconvolution reference matrix is deconvoluted for a second time to obtain a final cell type proportion and expression. Compared with the previous deconvolution method, the embodiments of the present specification obtain an initial cell type proportion by deconvoluting based on the initial reference matrix, solve the problem of incomplete cell types in the traditional reference matrix, and remove the noise in the initial reference matrix through the update of the reference matrix, thereby reducing the influence of certain redundant cell types in the initial reference matrix on the analysis result, and improving the accuracy of deconvolution.

[0058] The above description is only a summary of some technical schemes of the embodiments of the present specification, in order to more clearly understand the technical means of some embodiments of the present specification, and the contents of the specification can be implemented, and in order for the above and other purposes, features and advantages of the embodiments of the present specification to be more obvious and easy to understand, the following preferred embodiments are specifically described below, and the accompanying drawings are described in detail as follows. BRIEF DESCRIPTION OF DRAWINGS

[0059] In order to more clearly illustrate the technical schemes in the embodiments of the present specification or the prior art, the following will briefly introduce the drawings needed to be used in the embodiment or prior art description. Obviously, the drawings in the following description are only some embodiments of the present specification, and those skilled in the art can also obtain other drawings according to these drawings without creating any creative labor.

[0060] Figure 1 A step schematic diagram of a bulk data cell analysis method in some embodiments of the present specification is shown;

[0061] Figure 2 A block diagram of a bulk data cell analysis method in some embodiments of the present specification is shown;

[0062] Figure 3 A schematic diagram of an initial reference matrix in some embodiments of the present specification is shown;

[0063] Figure 4 A step schematic diagram of constructing an initial reference matrix in some embodiments of the present specification is shown;

[0064] Figure 5 Figure 8(a) shows a true value diagram of a data set in some embodiments of the present specification;

[0065] Figure 6 Figure 8(b) shows a diagram of a second deconvolution result in some embodiments of the present specification;

[0066] Figure 7 Figure 8(c) shows a diagram of a third deconvolution result in some embodiments of the present specification;

[0067] Figure 8(a) shows a true value diagram of a data set in some embodiments of the present specification;

[0068] Figure 8(b) shows a diagram of a second deconvolution result in some embodiments of the present specification;

[0069] Figure 9 Figure 8(c) shows a diagram of a third deconvolution result in some embodiments of the present specification;

[0070] Figure 10 Figure 8(a) shows a true value diagram of a data set in some embodiments of the present specification;

[0071] Figure 11 Figure 8(b) shows a diagram of a second deconvolution result in some embodiments of the present specification;

[0072] BRIEF DESCRIPTION OF THE DRAWINGS

[0073] 1001, an acquisition module;

[0074] 1002, a construction module;

[0075] 1003, a first deconvolution module;

[0076] 1004, an update module;

[0077] 1005, a second deconvolution module;

[0078] 1102, a computer device;

[0079] 1104, a processor;

[0080] 1106, a memory;

[0081] 1108, a driving mechanism;

[0082] 1110, an input / output module;

[0083] 1112, an input device;

[0084] 1114, an output device;

[0085] 1116, presentation device;

[0086] 1118, graphical user interface;

[0087] 1120, network interface;

[0088] 1122, communication link;

[0089] 1124, communication bus. DETAILED DESCRIPTION

[0090] The technical solutions in the embodiments of the present specification will be described clearly and completely below in combination with the drawings in the embodiments of the present specification. Obviously, the described embodiments are only some of the embodiments of the present specification, but not all the embodiments. Based on the embodiments in the present specification, all other embodiments obtained by those skilled in the art without creative work are within the scope of protection of the present specification.

[0091] It should be noted that the terms "first", "second", and the like in the present specification and claims and the above-described drawings are used to distinguish similar objects, and do not necessarily have to describe a specific order or a chronological order. It should be understood that the data used in this way can be exchanged under appropriate circumstances, so that the embodiments of the present specification described herein can be implemented in an order other than those illustrated or described herein. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, device, product or apparatus that includes a series of steps or units does not have to be limited to those steps or units clearly listed, but can include other steps or units that are not clearly listed or inherent to these processes, methods, products or apparatus.

[0092] It should be noted that the user information (including but not limited to user equipment information, user personal information, etc.) and data (including but not limited to data for analysis, stored data, displayed data, etc.) involved in the present application are all information and data authorized by the user or authorized by all parties. And the acquisition, storage, use, processing, etc. of data in the technical solutions described in the embodiments of the present application comply with relevant regulations.

[0093] In order to solve the above problems, the embodiments of the present specification provide a bulk data cell analysis method, Figure 1is a schematic diagram of a step of a bulk data cell analysis method provided by an embodiment of the present specification. The present specification provides method operation steps as described in the embodiments or flowcharts, but more or fewer operation steps can be included based on conventional or non-inventive labor. The order of steps listed in the embodiments is only one of the many step execution orders, and does not represent the only execution order. In actual system or device product execution, the method order shown in the embodiments or the drawings can be executed in sequence or in parallel. Specifically, as shown in Figure 1 and Figure 2 The method can include:

[0094] S101: Obtain a single-cell data set and a cell type annotation template, the cell type annotation template including a pre-defined cell type hierarchy and a feature gene set.

[0095] It can be understood that single-cell data refers to gene expression data at the level of individual cells obtained through single-cell sequencing technology. In the embodiments of the present specification, single-cell data is obtained from the GEO database. The cell type annotation template is a standard data set for cell type annotation of single-cell data, which defines 20 non-cancer cell types, cancer cell types and their hierarchical relationships. The hierarchical relationships can be represented by parent nodes and child nodes, such as immune cells as parent nodes, including T cells, B cells and other child nodes, and T cells can be further divided into CD4+ T and CD8+ T cell child nodes. The feature gene set is a list of characteristic genes specific to each cell type, such as CD4 gene for CD4+ T cells and CD8A gene for CD8+ T cells. In the embodiments of the present specification, the cell type annotation template can be obtained from the scATOMIC model.

[0096] S102: Construct an initial reference matrix according to the single-cell data set and the cell type annotation template, each element in the initial reference matrix representing the gene expression amount of each cell state under each feature gene.

[0097] In the embodiments of the present specification, each single cell can be labeled with its belonging cell type / state through the cell type annotation template, forming an annotated single-cell data set, wherein the cell state is a subdivision of the cell type. According to the annotated cell type / state in the single-cell data set, the expression value of the feature gene in each cell state can be calculated, thereby constructing the initial reference matrix, as shown in Figure 3The initial reference matrix is shown in FIG. 1. The rows of the initial reference matrix represent signature genes, the columns represent cell states, and the matrix elements represent the expression level of the signature genes in the corresponding cell state. When only one type of cancer cell is analyzed, the initial reference matrix contains 21 cell types, i.e. 20 non-cancer cell types and one specific cancer cell type. When multiple cancer analysis such as metastatic cancer analysis is performed, the initial reference matrix contains more than one cancer cell type and 20 non-cancer cell types. The 20 cell types include hematopoietic stem / progenitor cells, mast cells, B cells, plasma cells, monocytes, macrophages, NK cells, CD4+ T cells, CD8+ T cells, dendritic cell specific subtype (ASDC), conventional dendritic cells, plasmacytoid dendritic cells, oligodendrocytes, enteric glial cells, smooth muscle cells, endothelial cells, fibroblasts, myofibroblasts, cancer-associated fibroblasts, and cancer-associated myofibroblasts.

[0098] S103: deconvolving the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result, wherein the first deconvolution result comprises a first proportion and a first gene expression level of each cell type.

[0099] In the embodiments of the present specification, the bulk data to be analyzed refers to the gene expression data of a tissue or sample obtained by a conventional sequencing technology (such as RNA-seq), which reflects the overall gene expression level of a mixed cell population, wherein the heterogeneity of individual cells cannot be distinguished. In the embodiments of the present specification, the bulk data to be analyzed is one or more. By deconvolving the bulk data to be analyzed, the proportion of cell types contained in the bulk data and the gene expression level of each cell type can be inferred, i.e. the mixed signal is decomposed into single components, which is of great significance for analyzing the composition of immune cells in the tumor microenvironment and the cell heterogeneity in the tissue sample.

[0100] S104: updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix.

[0101] It can be understood that by using the information in the first deconvolution result, such as the inferred gene expression of cell types, the initial reference matrix is iteratively optimized to form a first reference matrix that is closer to the bulk data to be analyzed, thereby reducing the deviation of the initial reference matrix from the actual bulk data and improving the accuracy of subsequent deconvolution.

[0102] S105: deconvolving the bulk data to be analyzed according to the first reference matrix to obtain a second deconvolution result, wherein the second deconvolution result comprises a second proportion and a second gene expression level of each cell type.

[0103] It can be understood that using the updated first reference matrix to deconvolute the bulk data to be analyzed again can obtain more accurate cell type proportions and expression amounts. Through iterative optimization, the second deconvolution result is usually closer to the real cell composition and expression pattern than the first deconvolution result, and is especially suitable for scenarios where the reference matrix and the target data have batch effects or sample specificity.

[0104] With the technical solutions described above, the bulk data cell analysis method provided by the embodiments of the present specification first obtains an initial reference matrix containing multiple cell types through a single-cell data set and a cell type annotation template, the initial reference matrix is composed of a pan-cancer tumor microenvironment (TME) cell hierarchy, each parent node in the cell hierarchy represents a group of related cells, and each terminal node represents a target single-cell category. Then, the bulk data to be analyzed is deconvoluted based on a deconvolution algorithm and the initial reference matrix to obtain a first deconvolution result, and then an updated deconvolution reference matrix is obtained according to the first deconvolution result. The updated deconvolution reference matrix is deconvoluted a second time to obtain final cell type proportions and expressions. Compared with previous deconvolution methods, the embodiments of the present specification obtain initial cell type proportions based on the initial reference matrix, solve the problem of incomplete cell types in traditional reference matrices, and through updating of the reference matrix, remove noise in the initial reference matrix and reduce the influence of certain redundant cell types in the initial reference matrix on the analysis result, thereby improving the accuracy of deconvolution.

[0105] In the embodiments of the present specification, the initial reference matrix is constructed according to the single-cell data set and the cell type annotation template, including: Figure 4

[0106] S401: annotating each single cell in the single-cell data set using the cell type annotation template to obtain a cell state label of each single cell;

[0107] S402: summing expression data of all single cells of the same cell state to obtain expression counts of each cell state under each characteristic gene;

[0108] S403: normalizing the expression counts of each cell state under each characteristic gene to obtain gene expression amounts of each cell state under each characteristic gene;

[0109] S404: constructing an initial reference matrix according to the gene expression amounts of each cell state under each characteristic gene.

[0110] ​Specifically, the original expression amount of each single cell c in each gene g is obtained, the gene expression data of all single cells of the same cell state is summed based on the cell state, and the total expression count of the cell state can be obtained. For example, 3 cell states c1, c2 and c3 are defined in the annotation template, and each state contains 2 single cells: single cell 1 and single cell 2, and the characteristic gene set includes two genes: gene A and gene B, and the gene expression amount of each single cell under each cell state is known, then the expression count of each cell state under each gene and the total count can be calculated, for example, the total count of gene A of c1 is obtained by summing the counts of gene A of single cell 1 and single cell 2 in c1, and the total count of gene B of c1 is obtained in the same way, and thus the total expression count of c1 is obtained. The total count of gene A of c1 is normalized according to the total expression count, and the gene expression amount of c1 under gene A is obtained. Taking gene A and gene B as the rows of the matrix and the cell state as the column of the matrix, the gene expression amount of each cell state under each gene is taken as the matrix element, and the initial reference matrix is constructed In this way, by converting the gene expression information at the single cell level into the characteristic matrix of the cell state, prior knowledge can be provided for the subsequent deconvolution algorithm, thereby supporting the accurate inference of the cell composition in the bulk data.

[0111] In the embodiments of the present specification, with reference to Figure 5 , the deconvolution of the initial reference matrix on the bulk data to be analyzed is performed to obtain a first deconvolution result, comprising:

[0112] S1: using Gibbs sampling method to determine the joint posterior distribution of cells in the bulk data to be analyzed by taking the initial reference matrix and the bulk data to be analyzed as the prior distribution and sample information respectively;

[0113] S2: marginalizing the joint posterior distribution to obtain the gene expression marginal distribution matrix and the proportion marginal distribution of the cell state in the bulk data to be analyzed;

[0114] S3: calculating the gene expression mean matrix and the proportion mean of the cell state according to the gene expression marginal distribution matrix and the proportion marginal distribution of the cell state;

[0115] S4: summing the gene expression mean matrix and the proportion mean of all cell states belonging to the same cell type to obtain the gene expression matrix and the proportion of the cell type;

[0116] S5: updating the initial reference matrix according to the gene expression matrix of the cell type to obtain a fourth reference matrix;

[0117] S6: taking the fourth reference matrix as a new prior distribution, iteratively performing steps S1-S5 until a preset iteration termination condition is reached, to obtain a first deconvolution result.

[0118] It can be understood that, in the embodiments of the present specification, the BayesPrism deconvolution algorithm is used to infer the cell type proportion and expression of the bulk data to be analyzed based on the initial reference matrix as prior knowledge. The initial reference matrix and the expression X of the input nth bulk data to be analyzed are taken as prior distribution and sample information respectively, and the Gibbs sampling method is used to iteratively sample the input data, to infer the joint posterior distribution of the cell state proportion μ n and the gene expression U n of the bulk data. n For the nth bulk data, the marginal distribution of the gene expression matrix and the proportion of the cell state is obtained by marginalizing the joint posterior from the joint posterior distribution and , so as to estimate the gene expression matrix U n and the proportion μ n of each cell state, and the point estimate is obtained by averaging the sampling results of the marginal distribution and as the best estimate of the expression and proportion of the cell state.

[0119] If the same cell type contains multiple cell states, such as CD4+ T cells containing multiple subtypes, the posterior values of the cell states in each cell type and are summed up, and all cell states of each cell type are combined to estimate the gene expression matrix Z n and the proportion θ 0n of each cell type, that is θ 0n ​The initial cell type proportion of the nth bulk data. The cell types can be broadly divided into non-malignant cells (such as immune cells, stromal cells) and malignant cells (such as tumor cells), and the two types have large differences in gene expression patterns. Due to gene mutation, chromosomal abnormality, epigenetic remodeling, etc., the gene expression between different individuals, different lesions, and even cells in the same lesion is significantly different. In the bulk data of tumor tissue, malignant cells usually account for a dominant proportion, and their expression signals may mask the weak signals of non-malignant cells. Non-malignant cells belong to normal somatic cells and have relatively stable core marker gene expression. In the tumor microenvironment, the proportion of non-malignant cells is usually low, and their expression signals are easily masked by malignant cells. Directly inferring from the bulk data is susceptible to noise. If a fixed reference matrix is used for deconvolution, the heterogeneity between samples will be ignored, resulting in bias in proportion and expression estimation. Therefore, the reference matrix of malignant cells and non-malignant cells in each bulk data is inferred respectively, and the initial reference matrix is dynamically updated to make it more suitable for the cell characteristics of the current sample. In the embodiments of the present specification, the initial reference matrix is updated according to the gene expression matrix of the cell type to obtain a fourth reference matrix, including: Figure 6

[0120] S601: A second reference matrix of malignant cells and an initial reference matrix of non-malignant cells are calculated according to the gene expression matrix of the cell type.

[0121] S602: A third reference matrix of non-malignant cells is calculated according to the initial reference matrix of non-malignant cells and a preset prior probability.

[0122] S603: The second reference matrix and the third reference matrix are merged to obtain a fourth reference matrix.

[0123] Specifically, in step S601, first, the gene expression matrix Z of the cell type of the bulk data to be analyzed is obtained. n The expression value of the malignant cell under its characteristic gene is inferred by Gibbs sampling, and then the expression estimation value is taken as the element of the matrix, the characteristic gene is taken as the row, and the malignant cell type is taken as the column. The second reference matrix ψ is constructed. mal In the embodiments of the present specification, the expression estimation value of each characteristic gene of the malignant cell is calculated by the following formula:

[0124]

[0125] Wherein, represents the expression value of the characteristic gene g of the malignant cell, t1 represents the malignant cell, represents the expression estimation value of the gene g of the malignant cell in the gene expression matrix. ​

[0126] For the third reference matrix of non-malignant cells, in the embodiments of the present specification, the original counts of all cells of gene g in non-malignant cell type t2 are first summed, that is, the total expression amount of non-malignant cell type t2 at gene g, and then the original counts of all genes of all cells in non-malignant cell type t2 are summed to obtain the overall sequencing depth of non-malignant cell type t2, and the relative expression proportion of gene g in cell type t2 can be calculated according to the above two values, and the initial reference matrix of non-malignant cells can be constructed according to the relative expression proportion wherein the relative expression proportion of gene g in cell type t2 can be expressed by the following formula:

[0127]

[0128] wherein, represents the matrix value in the initial reference matrix of non-malignant cells, t2 represents non-malignant cells, W cg represents the original count of the cth cell at the gth gene. Since the content of non-malignant cells is relatively small, the present specification sets a prior probability conforming to a normal distribution as one of the calculation parameters, that is, to improve the estimation accuracy of low-abundance cells. In the embodiments of the present specification, the matrix value of each feature gene of non-malignant cells in the reference matrix is calculated according to the initial reference matrix of non-malignant cells and the preset prior probability by the following formula:

[0129]

[0130] wherein, represents the matrix value in the reference matrix of each feature gene of non-malignant cells, represents the matrix value in the initial reference matrix of non-malignant cells, represents the prior probability. The reference matrix of each feature gene of non-malignant cells can be constructed according to the matrix value in the reference matrix of each feature gene of non-malignant cells, and finally the reference matrix of each feature gene of non-malignant cells is merged to obtain the third reference matrix ψ of non-malignant cells env .

[0131] In the embodiments of the present specification, the second reference matrix ψ mal and the third reference matrix ψ env can be obtained by merging the updated reference matrix ψ n , and ψ n is taken as the latest prior distribution, that is, to re-estimate the cell type proportion and the marginal posterior of gene expression of each bulk sample, wherein is the cell type proportion of the updated reference matrix, until the parameter estimation converges, such as the cell proportion and expression amount no longer significantly changes.

[0132] In the embodiments of the present specification, with reference to Figure 7 , the updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix comprises:

[0133] S701: excluding cell types with a first proportion less than a first preset threshold in the first deconvolution result to obtain a first cell type set;

[0134] S702: excluding cell types with a correlation less than a second preset threshold with a cancer type in the first cell type set to obtain a second cell type set;

[0135] S703: obtaining a first reference matrix according to the first gene expression amount of each cell type in the second cell type set.

[0136] It can be understood that in the first deconvolution result, the cell types with small proportions and the cell types with small correlations with cancer are deleted in the initial reference matrix, the reference matrix is re-customized, and finally the updated first reference matrix is obtained. Based on the first reference matrix, the second deconvolution is performed to obtain the final result. In this way, low-abundance cell types and cell types unrelated to cancer can be filtered out, so as to reduce noise in the deconvolution process, focus on the main components, and improve the deconvolution accuracy. To verify the effectiveness of the bulk data analysis in the embodiments of the present specification, the pre-constructed pseudo bulk data is used as test data, and the deconvolution result of the embodiments of the present specification is evaluated by the root mean square error, the Pearson correlation coefficient and the cosine similarity. FIG. 8(a) is a ground truth diagram of the data set used, and FIG. 8(b) is the deconvolution result obtained by using the embodiments of the present specification. It can be seen that the deconvolution result obtained by the embodiments of the present specification has a small difference from the ground truth in the data set. The embodiments of the present specification are also compared with other deconvolution methods, and the comparison results of the model effects are shown in Figure 9 Figure 9 The deconvolution result diagrams of six existing deconvolution methods are included in Figure 9 ​The deconvolution result in the above table is compared with FIG. 8(b), and it can be found that the methods of CIBERSORT, xcell, eqic and quantiseq all have the problem of missing cell types, and the effect is poor compared with the initial reference matrix containing all cell types proposed in the embodiments of the present specification. Although the Bisque method uses the pan-cancer reference matrix of the embodiments of the present specification for deconvolution, it can include all cell types, but the overall effect is poor, because the model algorithm does not consider the problem of cell heterogeneity, and the embodiments of the present specification update the reference matrix and update the template in the algorithm, which removes the noise effect and improves the ability of the model to process cell heterogeneity.

[0137] Based on the bulk data cell analysis method described above, the embodiments of the present specification also correspondingly provide a bulk data cell analysis device. The device can include a system (including a distributed system), software (application), module, component, server, client, etc. using the method described in the embodiments of the present specification, and a device combined with necessary implementation hardware. Based on the same innovative concept, the device in one or more embodiments provided by the embodiments of the present specification is as described in the following embodiments. Since the implementation scheme of the device to solve the problem is similar to the method, the implementation of the specific device of the embodiments of the present specification can be referred to the implementation of the foregoing method, and the repeated parts will not be described herein. The term "unit" or "module" used below can be a combination of software and / or hardware that implements a predetermined function. Although the device described in the following embodiments is preferably implemented in software, hardware, or a combination of software and hardware is also possible and is conceived.

[0138] Specifically, Figure 10 is a schematic diagram of the module structure of one embodiment of the bulk data cell analysis device provided by the embodiments of the present specification, referring to Figure 10 The bulk data cell analysis device provided by the embodiments of the present specification includes:

[0139] The acquisition module 1001 is configured to acquire a single-cell data set and a cell type annotation template, wherein the cell type annotation template includes a pre-defined cell type hierarchy and a feature gene set;

[0140] The construction module 1002 is configured to construct an initial reference matrix according to the single-cell data set and the cell type annotation template, wherein each element in the initial reference matrix represents a gene expression amount of each cell state under each feature gene;

[0141] The first deconvolution module 1003 is configured to perform deconvolution on the bulk data to be analyzed according to the initial reference matrix, to obtain a first deconvolution result, wherein the first deconvolution result includes a first proportion and a first gene expression amount of each cell type;

[0142] The updating module 1004 is configured to update the initial reference matrix according to the first deconvolution result to obtain a first reference matrix.

[0143] The second deconvolution module 1005 is configured to perform deconvolution on the bulk data to be analyzed according to the first reference matrix to obtain a second proportion and a second gene expression amount of each cell type in the bulk data to be analyzed.

[0144] The device provided by the embodiments of the present specification has the same beneficial effects as the method provided by the present specification, which will not be repeated here.

[0145] Referring to Figure 11 As shown in FIG. 11, based on the bulk data cell analysis method described above, an embodiment of the present specification further provides a computer device 1102, wherein the method described above runs on the computer device 1102. The computer device 1102 can include one or more processors 1104, such as one or more central processing units (CPUs), each of which can implement one or more hardware threads. The computer device 1102 can also include any memory 1106 for storing any kind of information, such as code, settings, data, etc. Without limitation, for example, the memory 1106 can include any one or a combination of more than one of the following: any type of RAM, any type of ROM, a flash memory device, a hard disk, an optical disk, etc. More generally, any memory can use any technology for storing information. Further, any memory can provide volatile or non-volatile retention of information. Further, any memory can represent a fixed or removable component of the computer device 1102. In one case, the computer device 1102 can perform any operation of the associated instructions when the processor 1104 executes the associated instructions stored in any memory or combination of memories. The computer device 1102 also includes one or more drive mechanisms 1108, such as a hard disk drive mechanism, an optical disk drive mechanism, etc., for interacting with any memory.

[0146] The computer device 1102 can also include an input / output module 1110 (I / O) for receiving input (via input device(s) 1112) and for providing output (via output device(s) 1114). One specific output mechanism can include a presentation device 1116 and associated graphical user interface (GUI) 1118. In other embodiments, the input / output module 1110 (I / O), the input device(s) 1112, and the output device(s) 1114 can not be included, e.g., as merely a computer device in a networked environment. The computer device 1102 can also include one or more network interfaces 1120 for exchanging data with other devices via one or more communication links 1122. The one or more communication buses 1124 couple the various components in the computer device 1102.

[0147] The communication links 1122 can be any communication links, such as local area networks, wide area networks, e.g., the Internet, point-to-point connections, etc., or any combination thereof. The communication links 1122 can include any combination of hardwired links, wireless links, routers, gateway functionality, name servers, etc. governed by any protocol or combination of protocols.

[0148] Corresponding to the method as shown in Figures 1 to 7 The embodiment of the present specification also provides a computer readable storage medium, and the computer readable storage medium stores a computer program. The computer program is run by a processor to execute the steps of the method.

[0149] The embodiment of the present specification also provides a computer readable instruction, wherein the program in the computer readable instruction makes the processor execute the method as shown in Figures 1 to 7 .

[0150] The embodiment of the present specification also provides a computer program product, including at least one instruction or at least one program. The at least one instruction or the at least one program is loaded and executed by a processor to implement the method as shown in Figures 1 to 7 .

[0151] It should be understood that the size of the sequence number of each process described above in various embodiments of the present specification does not mean the order of execution, and the execution order of each process should be determined according to its function and inherent logic, and should not constitute any limitation on the implementation process of the embodiment of the present specification.

[0152] It should also be understood that in the embodiment of the present specification, the term "and / or" is only a description of the association relationship of the associated objects, which means that there can be three relationships. For example, A and / or B can represent the following three cases: A exists alone, A and B exist together, and B exists alone. In addition, the character " / " in the present specification generally represents an "or" relationship between the front and rear associated objects.

[0153] Those skilled in the art can clearly understand that the units and algorithm steps of each example described in combination with the embodiments disclosed in the specification can be realized in electronic hardware, computer software or a combination of both. In order to clearly illustrate the interchangeability of hardware and software, the components and steps of each example have been described in the foregoing description in a general manner. Whether the functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. A person skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the specification.

[0154] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working processes of the above-described system, device and unit can refer to the corresponding processes in the foregoing method embodiments, which will not be repeated here.

[0155] In several embodiments provided in the specification, it should be understood that the disclosed system, device and method can be implemented in other ways. For example, the above-described device embodiments are only schematic, for example, the division of the units is only a logical function division, and actual implementation can have another division manner, for example, a plurality of units or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the units shown or discussed can be indirect coupling or communication connection through some interface, device or unit, and can also be electrical, mechanical or other forms of connection.

[0156] The units described as separate components can or can not be physically separated, and the components shown as units can or can not be physical units, that is, they can be located in one place, or can be distributed on a plurality of network units. Part or all of the units can be selected according to actual needs to achieve the purpose of the embodiment of the specification.

[0157] In addition, each functional unit in each embodiment of the specification can be integrated in one processing unit, or each unit can be physically present separately, or two or more units can be integrated in one unit. The integrated unit can be realized in the form of hardware or in the form of a software functional unit.

[0158] The integrated unit, if implemented in the form of a software function unit and sold or used as an independent product, can be stored in a computer readable storage medium. Based on such understanding, the technical solutions of the specification or the entire or part of the technical solutions that essentially contribute to the prior art can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a plurality of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in the various embodiments of the specification. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), a magnetic disk or an optical disk, and various media that can store program codes.

[0159] The principles and implementation manners of the specification are described by using specific embodiments in the specification. The above embodiment description is only used to help understand the method of the specification and its core idea; meanwhile, for those skilled in the art, according to the idea of the specification, the specific implementation manners and application ranges will have changes. In summary, the content of the specification should not be understood as a limitation of the specification.

Claims

1. A cell analysis method for bulk data, characterized in that: The method comprises: Obtaining a single-cell dataset and a cell type annotation template, wherein the cell type annotation template includes a predefined cell type hierarchy and a set of characteristic genes; Constructing an initial reference matrix based on the single-cell dataset and the cell type annotation template, wherein each element in the initial reference matrix represents the gene expression level of each cell state under each characteristic gene; Deconvolving the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result, wherein the first deconvolution result includes a first ratio and a first gene expression level of each cell type; Updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix; Deconvolution is performed on the bulk data to be analyzed according to the first reference matrix to obtain a second deconvolution result, where the second deconvolution result includes a second proportion and a second gene expression level of each cell type.

2. The method according to claim 1, characterized in that The constructing of an initial reference matrix based on the single cell dataset and the cell type annotation template comprises: Annotating each single cell in the single cell dataset using the cell type annotation template to obtain a cell state label for each single cell; According to the cell state label of each single cell, the expression data of all single cells in the same cell state are summed to obtain the expression count of each cell state under each characteristic gene; Normalize the expression counts of each cell state under each characteristic gene to obtain the gene expression value of each cell state under each characteristic gene; An initial reference matrix was constructed based on the gene expression levels of each cell state under each characteristic gene.

3. The method according to claim 1, characterized in that The step of performing deconvolution on the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result includes: S1: using the initial reference matrix and the bulk data to be analyzed as the prior distribution and sample information, respectively, and using the Gibbs sampling method to determine the joint posterior distribution of cells in the bulk data to be analyzed; S2: marginalizing the joint posterior distribution to obtain the gene expression marginal distribution matrix and proportional marginal distribution of the cell state in the bulk data to be analyzed; S3: calculating the gene expression mean matrix and proportional mean of the cell state according to the gene expression marginal distribution matrix and proportional marginal distribution of the cell state; S4: Sum the gene expression mean matrix and ratio mean of all cell states belonging to the same cell type to obtain the gene expression matrix and ratio of the cell type; S5: updating the initial reference matrix according to the gene expression matrix of the cell type to obtain a fourth reference matrix; S6: Using the fourth reference matrix as a new prior distribution, iterate steps S1-S5 until a preset iteration termination condition is reached to obtain a first deconvolution result.

4. The method according to claim 3, characterized in that The cell types include malignant cells and non-malignant cells, and the initial reference matrix is ​​updated according to the gene expression matrix of the cell type to obtain a fourth reference matrix, including: According to the gene expression matrix of the cell type, a second reference matrix of malignant cells and an initial reference matrix of non-malignant cells are calculated; calculating a third reference matrix of non-malignant cells based on the initial reference matrix of non-malignant cells and a preset prior probability; The second reference matrix and the third reference matrix are combined to obtain a fourth reference matrix.

5. The method according to claim 4, characterized in that The step of calculating a second reference matrix of malignant cells based on the gene expression matrix of the cell type comprises: According to the gene expression matrix of the cell type, the expression value of each characteristic gene of the malignant cell is calculated using the following formula: in, Represents the second reference matrix of malignant cells in the bulk data to be analyzed, t1 represents malignant cells, Represents the gene expression matrix of malignant cells in the bulk data to be analyzed; A second reference matrix is ​​constructed according to the expression value of each characteristic gene of the malignant cells.

6. The method according to claim 4, characterized in that Based on the gene expression matrix of the cell type, an initial reference matrix of non-malignant cells is calculated, including: Based on the gene expression matrix of the cell type, the matrix value in the initial reference matrix of non-malignant cells is calculated using the following formula: in, represents the matrix value in the initial reference matrix of non-malignant cells, t2 represents non-malignant cells, W cg represents the raw count of the cth cell under the gth gene in non-malignant cells; An initial reference matrix is ​​constructed according to the matrix values ​​in the initial reference matrix.

7. The method according to claim 4, characterized in that The step of calculating a third reference matrix of non-malignant cells based on the initial reference matrix of non-malignant cells and a preset prior probability includes: The matrix value in the reference matrix of each characteristic gene of non-malignant cells was calculated using the following formula: in, Represents the matrix value in the reference matrix of each characteristic gene of non-malignant cells, represents the matrix values ​​in the initial reference matrix of non-malignant cells, represents the prior probability; Constructing a reference matrix for each characteristic gene of non-malignant cells according to the matrix values ​​in the reference matrix for each characteristic gene of non-malignant cells; The reference matrices of each characteristic gene of non-malignant cells are merged to obtain a third reference matrix of non-malignant cells.

8. The method according to claim 1, characterized in that The updating the initial reference matrix according to the first deconvolution result to obtain a first reference matrix includes: Screening out cell types whose first proportion in the first deconvolution result is less than a first preset threshold to obtain a first cell type set; Eliminating cell types whose correlation with cancer types is less than a second preset threshold from the first cell type set to obtain a second cell type set; A first reference matrix is ​​obtained according to the first gene expression level of each cell type in the second cell type set.

9. A bulk data cell analysis device, characterized in that: The device comprises: An acquisition module is used to acquire a single-cell dataset and a cell type annotation template, wherein the cell type annotation template includes a predefined cell type hierarchy and a set of characteristic genes; A construction module is used to construct an initial reference matrix based on the single-cell dataset and the cell type annotation template, wherein each element in the initial reference matrix represents the gene expression level of each cell state under each characteristic gene; a first deconvolution module, configured to perform deconvolution on the bulk data to be analyzed according to the initial reference matrix to obtain a first deconvolution result, wherein the first deconvolution result includes a first ratio and a first gene expression level of each cell type; An updating module, configured to update the initial reference matrix according to the first deconvolution result to obtain a first reference matrix; The second deconvolution module is configured to perform deconvolution on the bulk data to be analyzed according to the first reference matrix to obtain a second deconvolution result, where the second deconvolution result includes a second proportion and a second gene expression level of each cell type.

10. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the method according to any one of claims 1 to 8 is implemented.

11. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the method according to any one of claims 1 to 8 is implemented.

12. A computer program product, characterized in that The method comprises at least one instruction or at least one program, wherein the at least one instruction or the at least one program is loaded and executed by a processor to implement the method according to any one of claims 1 to 8.

Images