Cell fate mapping research method based on CytoEnergy technology

Abstract

The present invention provides a cell fate mapping research method based on CytoEnergy technology, which belongs to the field of biorecognition technology, and includes the following steps: step S1, preprocessing single cell transcriptome data; step S2, constructing a local gene network based on the Hopfield neural network; step S3, calculating the local gene network energy and cell energy, and normalizing the cell energy; step S4, introducing a Gaussian mixture model, identifying the potential state of the cell by analyzing the energy distribution difference of the cell, and clustering cells in the same differentiation potential state; step S5, identifying cells in different fate terminal states; step S6, inferring the cell lineage trajectory. Through the above method, without relying on prior knowledge, the technology provides directional information for de novo inference of cell fate transformation from the perspective of cell plasticity, and can accurately restore the dynamic change process of cell fate.

Description

Technical Field

[0001] The present invention relates to the field of biorecognition technology, and in particular to a cell fate mapping research method based on CytoEnergy technology. Background Art

[0002] Cell fate mapping describes the transition process of cells from one state to another, and is particularly critical for revealing how cells respond to changes in the internal and external environment. In multicellular organisms, the formation of complex tissue structures depends on various types of cells, which show essential differences in function and plasticity. These cells undergo a series of state transitions during various processes such as development, differentiation, and reprogramming, which usually occur in a highly asynchronous manner and are accompanied by changes in epigenetic and transcriptional levels. Although advances in single-cell sequencing technology have made it possible to capture these dynamic changes, the destructiveness of cell lysis during the sequencing process makes it impossible to track the dynamic changes of cell states in real time, and only a series of static snapshots of cell fate transitions can be obtained.

[0003] To overcome this problem, a variety of algorithmic frameworks have been developed to map the dynamic characteristics of cell fate from discrete sequencing data. However, most of the existing trajectory inference methods (such as Monocle 1 / 2, SCUBA, and Wanderlust) rely on prior biological knowledge (such as time information or specifying initial cells / states) to determine the direction of cell fate transformation. Lineage and clone tracking technologies provide important insights into the lineage fate of cells, but their wide application in vivo is limited due to the influence of cell type-specific genetic markers. In addition, methods based on RNA velocity estimation (such as Velocyto, scVelo, and Dynamo) infer the direction of cell fate by calculating the ratio of spliced ​​and unspliced ​​mRNA, but the noise introduced by velocity estimation and the inability to accurately capture the main cell fate trajectory are its biggest limitations. In this context, it is particularly critical to determine the direction of cell fate transformation through differences in cell plasticity. Although methods based on transcriptional diversity (CytoTRACE) and entropy (such as StemID, SLICE, and SCENT) provide directional information for cell fate transition by quantifying cell plasticity, some of these methods rely on additional external input information (such as protein-protein interaction network information), and some can only infer the direction of fate transition between different cell clusters, while ignoring the heterogeneity information within the cell clusters. Therefore, it is urgent to develop new technologies that can effectively restore the direction of cell fate transition and comprehensively map the dynamic changes of cell fate. Summary of the invention

[0004] The purpose of the present invention is to provide a cell fate mapping research method based on CytoEnergy technology. Without relying on prior knowledge, this technology provides directional information for de novo inference of cell fate transformation from the perspective of cell plasticity, and can accurately restore the dynamic change process of cell fate.

[0005] To achieve the above object, the present invention provides a cell fate mapping research method based on CytoEnergy technology, comprising the following steps:

[0006] Step S1, preprocessing single-cell transcriptome data;

[0007] Step S2, constructing a local gene network based on the Hopfield neural network;

[0008] Step S3, calculating the gene local network energy and cell energy, and normalizing the cell energy;

[0009] Step S4, introducing a Gaussian mixture model to identify the potential state of cells by analyzing the difference in cell energy distribution, and clustering cells in the same differentiation potential state;

[0010] Step S5, identifying cells in different terminal fate states;

[0011] Step S6: infer cell lineage trajectories.

[0012] Preferably, in step S1, the single cell transcriptome data is preprocessed, and the specific operations are:

[0013] Given a system containing N cells and M 1 Single-cell expression matrix of genes First, the single-cell transcriptome data was log-transformed using the NormalizeData function, and then standardized using the ScaleData function. Then, the highly variable genes were screened using the feature selection method. Finally, the highly variable genes were identified using the FindVariableFeatures function.

[0014] Processed dataset contains N cells and M highly variable genes, and the state of cell n is represented by the vector express; Represents the state value of gene i in cell n. The specific calculation formula is as follows:

[0015]

[0016] in, Represents the normalized expression level of gene i in cell n.

[0017] Preferably, in step S2, the nodes in the gene local network represent genes, and the interactions between genes in the network are specified by their adjacency matrix W, where the matrix elements w ij Represents the relationship between genes i and j:

[0018]

[0019] in, represents the state value of gene j in cell n; p represents the number of stable states.

[0020] Preferably, in step S3, the calculation formula of the gene local network energy is as follows:

[0021]

[0022] in, is the energy value of the local network centered on gene i in cell n; N(i) represents the neighborhood of node i in the network.

[0023] Preferably, in step S3, the cell energy consists of the local network energy levels of all genes contained in the cell, and is defined as:

[0024]

[0025] Among them, E n is the energy value of cell n; M represents the total number of highly variable genes.

[0026] Preferably, in step S3, the normalization process is:

[0027]

[0028] in, and Respectively represent the minimum and maximum energy values ​​of all cells; E' n is the normalized cellular energy.

[0029] Preferably, in step S4, a Gaussian mixture model is introduced to identify the potential state of cells by analyzing the energy distribution difference of cells, and cells in the same differentiation potential state are clustered. The specific operation is:

[0030] For a cell population S containing N cells, the normalized cell energy E' n∈N To quantify the differentiation potential status of each cell;

[0031] The cells are sorted according to their energy and the E' of all cells is calculated using a Gaussian mixture model. n∈N The values ​​are fitted;

[0032] Among them, the Gaussian mixture model is as follows:

[0033]

[0034] in, represents the kth Gaussian distribution with a given mean μ k and variance The probability density function under the condition λ k >0 is the mixing weight of the k-th Gaussian distribution, satisfying

[0035] The number of components K of the Gaussian distribution is determined by the Bayesian Information Criterion BIC, which is defined as:

[0036]

[0037] Here, v is the number of estimated parameters.

[0038] Preferably, in step S5, cells in different terminal fate states are identified by:

[0039] First, cell populations with different differentiation potentials were identified;

[0040] Next, all cells were ranked according to their relative positions in the entire data set;

[0041] When cells are arranged in chronological order from low to high in the actual biological process, the CytoEnergy technology defines the relative position ranking score of cells in the process of decreasing plasticity by Z'=1-E', where E' represents the normalized energy level of the cell; conversely, in the biological process of increasing cell plasticity, the position ranking of cells is determined by Z'=E';

[0042] CytoEnergy introduces a Gaussian distribution model, whose probability density function is defined as follows:

[0043]

[0044] Where μ represents Gaussian distribution The expected value of ; σ is the standard deviation, respectively expressed as:

[0045] μ=median(Z') (9)

[0046] σ=median(|Z'-μ|) (10)

[0047] Calculate Gaussian distribution The theta quantile of By default: θ = 85% where the position ranking score Z' exceeds the threshold cells, i.e. is considered to be in a potential terminal fate state; in biological processes with reduced plasticity, a low-potential state and Cells that are in a terminal fate are considered to be in a terminal fate state; in contrast, in biological processes with increased plasticity, cells that simultaneously satisfy a high-potential state and Cells with a terminal fate are considered to be in a terminal fate state, and the pseudotime value of cells in a terminal fate state is set to 1.

[0048] Preferably, in step S6, the cell lineage trajectory is inferred, and the specific operations are as follows:

[0049] The cell fate of the remaining cells is measured by the combined similarity index CSI, which is defined as follows:

[0050] CSI=(1-GD lt )*PC lt (11)

[0051] Among them, GD lt represents the Gower distance between the remaining cell l and the terminal state t, and is calculated as follows:

[0052]

[0053] in, represents the normalized expression level of gene i in cell l; is the average normalized expression level of gene i in all cells in terminal state t; d refers to the total number of genes in cell l;

[0054] PC lt represents the Pearson correlation coefficient between the remaining cells l and the terminal state t, which is defined as follows:

[0055]

[0056] According to the CSI maximization principle, each remaining cell l is assigned to the differentiation path P represented by the corresponding terminal state t t As shown below:

[0057]

[0058] Where T represents the total number of terminal cell fate states; L refers to the total number of remaining cells, c l Represents the remaining cells l; c t represents the terminal state t; the pseudo-time value of cell l is determined by its position score Z', which reflects the relative position of the cell in the fate trajectory;

[0059] The local weighted scatter point smoothing regression method was used, and the span parameter was set to 0.75 to fit the expression trend of genes along the lineage path;

[0060] A span parameter of 0.95 was set to analyze the changing trend of the local network energy of genes along the lineage trajectory.

[0061] Therefore, the present invention adopts the above-mentioned cell fate mapping research method based on CytoEnergy technology. CytoEnergy technology projects high-dimensional genetic data into a two-dimensional energy landscape to more intuitively display the energy state of cells. Without relying on prior knowledge, this technology provides directional information for de novo inference of cell fate changes from the perspective of cell plasticity, and can accurately restore the dynamic change process of cell fate. That is, CytoEnergy technology realizes unbiased and accurate cell fate dynamic mapping. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is a flow chart of a cell fate mapping research method based on CytoEnergy technology of the present invention;

[0063] Figure 2 This is an overview of pancreatic cell development on embryonic day 15.5 (E15.5) of mice; Figure 2 (A) is a diagram showing the relationship between the lineage development of pancreatic cells in the early embryo; Figure 2 (B) is a UMAP dimensionality reduction image showing the GLNE map of mouse pancreatic cell development E15.5 cells; Figure 2 (C) is a violin plot of the normalized energy values ​​of cells during mouse pancreatic cell development; Figure 2 (D) in the figure shows the energy density distribution during pancreatic development;

[0064] Figure 3 This is a reconstruction of the developmental pathway of pancreatic endocrine cells on day 15.5 of the mouse embryo; Figure 3 (A) shows the identification of terminal cells during pancreatic cell development; Figure 3 (B) in the figure shows the proportion of cells in different terminal states; Figure 3 (C) in the figure shows the identification of lineage pathways of different pancreatic endocrine cells; Figure 3 (D) in the middle is the pseudo-time trajectory analysis of pancreatic cell development;

[0065] Figure 4 is the dynamic pattern of key genes in the Alpha, Beta, Delta and Epsilon pathways; Figure 4 (A1) in the figure shows the GLNE level of the marker gene Gcg of Alpha cells in the four terminal cells; Figure 4(A2) shows the GLNE level of the Beta cell marker gene Ins2 in the four terminal cells; Figure 4 (A3) in the figure shows the GLNE level of Delta cell marker gene Sst in the four terminal cells; Figure 4 (A4) in the figure shows the GLNE level of the marker gene Ghrl of Epsilon cells in the four terminal cells; Figure 4 (B1) in the figure shows the changing trend of GLNE value of Gcg gene along four pathways; Figure 4 (B2) in the figure shows the changing trend of the GLNE value of the Ins2 gene along the four paths; Figure 4 (B3) in the figure shows the changing trend of the GLNE value of the Sst gene along the four paths; Figure 4 (B4) in the figure shows the changing trend of GLNE value of Ghrl gene along four pathways; Figure 4 (C1) in the figure shows the changing trend of the expression level of Gcg gene along the four pathways; Figure 4 (C2) in the figure shows the changing trend of the expression level of Ins2 gene along the four pathways; Figure 4 (C3) in the figure shows the changing trend of the expression level of Sst gene along the four pathways; Figure 4 (C4) in the figure shows the changing trend of the expression level of Ghrl gene along the four pathways. DETAILED DESCRIPTION

[0066] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0067] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0068] Embodiment 1

[0069] like Figure 1 As shown, the present invention provides a cell fate mapping research method based on CytoEnergy technology, including six parts: data preprocessing, quantification of gene local network energy, calculation of cell energy, identification of different potential states, recognition of terminal fate states, and inference of cell lineage trajectory.

[0070] Data preprocessing

[0071] In this example, single-cell transcriptome data of mouse pancreatic embryonic development E15.5 is involved. This dataset covers the dynamic development process of mouse pancreas in detail. Specifically, the unique molecular identifier count matrix comes from the dataset GSE132188.

[0072] First, the data were processed according to the preprocessing method described in the article from which the data came (Bastidas-Ponce A, Tritschler S, Dony L et al. Development. 2019 Jun 17; 146(12): dev173849). The final data obtained included 27,998 genes and 2,780 cells. These cells were annotated into different subpopulations ( Figure 2 (A)), including endocrine progenitor cells with low Ngn3 expression (Ngn3 low EP, n = 262), endocrine progenitor cells with high Ngn3 expression (Ngn3 high EP, n = 642), and Fev + Endocrine intermediate cells (Fev + , n = 592), as well as endocrine cells that produce different hormones, namely Alpha cells (glucagon-secreting cells; n = 481), Beta cells (insulin-secreting cells; n = 591), Delta cells (somatostatin-secreting cells; n = 70) and Epsilon cells (growth hormone-secreting cells; n = 142). Next, the above data were screened using feature selection methods. The specific operations are as follows: the data were logarithmically transformed and standardized using the NormalizeData and ScaleData functions of the Seurat software package, respectively. Finally, the highly variable genes were identified using the FindVariableFeatures function, where the parameter selection.method = "mvp". The state of cell n is represented by a vector express. Represents the status value of gene i in cell n, based on its expression level after standardized processing in cell n The sign of is determined as follows:

[0073]

[0074] Quantification of gene local network energy;

[0075] The gene local network was constructed with reference to the Hopfield Neural Network (HNN), which is a single-layer fully connected feedback neural network consisting of nodes and weighted undirected edges. The nodes of the network represent genes, and the interactions between genes are expressed through w ij To describe:

[0076]

[0077] in, and Represent the states of genes i and j in cell n respectively. p represents the number of stable states. Cells with clear phenotypic information are directly modeled as attractors (stable states) of HNN. For any given gene i in cell n, a local network with gene i as the core is constructed according to the method described in the previous section. The network consists of gene i and its first-order neighbor genes. The local network is given a local network energy The specific calculation formula is as follows:

[0078]

[0079] in, represents the state value of gene j in cell n, and N(i) represents the neighborhood of node i in the network (i.e., other nodes excluding node i). After traversing each gene, a GLNE matrix is ​​obtained. The elements of the matrix It represents the energy value of the local network centered on gene i in cell n. This value reflects the activity of gene i in cell n, where a higher GLNE value indicates a higher activity of the gene.

[0080] It is worth noting that the GLNE matrix E has the same dimension as the input gene expression matrix A, so this matrix was used for dimensionality reduction analysis in subsequent studies. The UMAP dimensionality reduction plot based on matrix E showed that mouse pancreatic endocrine cells exhibited four major fate lineage commitments during development: Alpha, Beta, Delta, and Epsilon, in which cells of the same type clustered together and followed similar energy characteristics ( Figure 2 (B) in the figure.

[0081] Calculation of cellular energy;

[0082] Energy value E of cell n n It consists of the local network energy levels of all genes it contains and is defined as:

[0083]

[0084] To facilitate subsequent comparative analysis, this embodiment normalizes the energy value of cell n, and the processed energy value E' n The definition is as follows:

[0085]

[0086] in, and Respectively represent the minimum and maximum energy values ​​of all cells. E' nThe value range is between 0 and 1. This value is used as a quantitative indicator of cell plasticity, reflecting the overall activity of the gene interaction network within each cell. A higher energy value reflects that the gene interaction network within the cell has a higher activity; conversely, the activity of the gene interaction network within the cell is lower.

[0087] According to the above steps, the energy of cells at different stages of pancreatic endocrine production in mice was calculated. The results showed that Ngn3 low EPs obtained the highest energy, followed by Ngn3 high EPs and Fev + In the intermediate stage, the energy level of the cell gradually decreases until it transforms into four types of endocrine cells, at which point the energy level of the cell reaches its lowest level ( Figure 2 Overall, the gradual decrease in energy levels during this process reflects the stepwise induction of endocrine cell fate.

[0088] Identification of different potential states;

[0089] CytoEnergy technology introduces a Gaussian Mixture Module (GMM) to identify the possible potential states of cells by analyzing the differences in cell energy distribution and clustering cells in the same potential state into one category. First, for a cell population S containing N cells, the normalized energy value E' n∈N Represents the differentiation potential of each cell. The energy value ranges from 0 to 1, reflecting the relative position of the cell in terms of potential state. Next, the cells are sorted according to the size of the energy value, and the E' of these cells is calculated using GMM. n∈N The value is fitted. This step aims to identify the different potential states existing in the cell population S and classify the cells according to the level of differentiation potential. The mathematical description of GMM is as follows:

[0090]

[0091] in, Indicates that the mean μ k and variance The probability density function of the k-th Gaussian distribution under λ. k >0 is the mixing weight of the k-th Gaussian distribution, satisfying The number of components K of the Gaussian distribution is determined by the Bayesian Information Criterion BIC, which is defined as follows:

[0092]

[0093] Where v is the number of estimated parameters, this step was implemented using the R package mclust (version 6.0.1). Specifically, the optimal number of cell potential states was determined by selecting the K value that maximizes the BIC.

[0094] Generally speaking, the potential states of cells can be divided into three categories: Low Potential State (LPS), the energy values ​​of these cells are at a lower level in the Gaussian distribution; High Potential State (HPS), the energy values ​​of these cells are at a higher level in the Gaussian distribution; Intermediate Potential State (IPS), the energy of these cells is between the two.

[0095] identification of terminal fate states;

[0096] In the absence of prior information, CytoEnergy technology can automatically identify cells in terminal fate states from single-cell gene expression data. The basic assumption of this technology is that in biological processes where cell plasticity decreases (such as cell differentiation), cells in terminal states usually show the lowest differentiation potential and are generally ranked higher than other cells in the entire trajectory. Conversely, in biological processes where cell plasticity increases (such as reprogramming), terminal state cells show the highest differentiation potential. In this case, cells with high differentiation potential (i.e., high energy) are selected as candidate cells in terminal fate states.

[0097] In order to effectively identify cells in a terminal fate state in a biological process where cell plasticity decreases, this embodiment first identifies a cell population with the lowest differentiation potential. This cell population is considered to be in a candidate terminal fate state and may be composed of multiple cell types. This step is performed according to the previous section. Next, all cells are sorted according to their relative position in the entire data set. Given that cells are arranged in chronological order from low to high in actual biological processes, CytoEnergy technology defines the relative position ranking score of cells in the process of decreasing plasticity by Z'=1-E', where E' represents the normalized energy level of the cell. Conversely, in a biological process where cell plasticity increases, the position ranking of cells is determined by Z'=E'.

[0098] Furthermore, CytoEnergy introduces a Gaussian distribution model, whose probability density function is defined as follows:

[0099]

[0100] Where μ represents Gaussian distribution The expected value of , σ is the standard deviation. They are expressed as:

[0101] μ=median(Z') (9)

[0102] σ=median(|Z'-μ|) (10)

[0103] The Gaussian distribution was further calculated The theta quantile of (by default θ = 85%), where the position ranking score Z' exceeds the threshold cells, i.e. It is considered to be in a potential terminal fate state. If both of the above conditions are met (i.e., in a specific potential state and the position ranking score exceeds ) are finally determined as terminal state cells, and the pseudo-time values ​​of these cells are set to 1.

[0104] According to the above steps, CytoEnergy technology successfully identified four endocrine terminal cell states, corresponding to Alpha, Beta, Delta and Epsilon cell states ( Figure 3 (A) in which Alpha and Beta cells occupy most of the terminal states, reflecting the main cell fates in pancreatic endocrine development ( Figure 3 (B) in the figure), a result that is consistent with known biological events. Due to the rarity of Delta cells in this data (3% of the total number in this data), it is difficult for RNA velocity-based methods to capture the differentiation state of Delta cells. For example, CellRank requires a stable macroscopic recovery of the Delta cell terminal state. scVelo is unable to effectively capture genes associated with Delta cell development. It is worth emphasizing that the CytoEnergy technology successfully captured the four main fate terminal states of pancreatic endocrine occurrence, of which 6.3% of the terminal cells were identified as Delta cells ( Figure 3 In summary, the CytoEnergy technology can effectively identify different terminal cell states, including Delta terminal cells, regardless of their proportion in the data.

[0105] Inference of cell lineage trajectories;

[0106] The previous steps successfully obtained a collection of cells in different terminal states {c (t) |t∈{1,2,…,T}}, where c (t) Represents cell fate t, which is composed of multiple cells in the same potential state; T represents the total number of terminal cell fate states. Next, we study how to assign the remaining cells to their respective cell fates. These remaining cells do not belong to the terminal state and are represented by l.

[0107] In order to solve the above problems, the combined similarity index CSI is introduced. CSI combines the correlation and the adjusted distance parameter to measure which different cell fates the remaining cells (cells that do not belong to the terminal state) are more inclined to transform into. The definition is as follows:

[0108] CSI=(1-GD lt )*PC lt (11)

[0109] Among them, GD lt represents the Gower distance between the remaining cell l and the terminal state t, and is calculated as follows:

[0110]

[0111] in, represents the normalized expression level of gene i in cell l, is the average normalized expression level of gene i across all cells in terminal state t, and the Pearson correlation coefficient between cell l and terminal state t is defined as follows:

[0112]

[0113] Next, according to the CSI maximization principle, each cell l is assigned to the differentiation path P represented by the corresponding terminal state t. t As shown below:

[0114]

[0115] During this process, P t It indicates that the cell starts from a high potential state (corresponding to a high energy level), passes through a series of transition cells in the intermediate state, and finally reaches an orderly path of different terminal states. These terminal states represent multiple possible fates of the cell. The pseudo-time value of cell l is determined by its position score Z', which reflects the relative position of the cell in the fate trajectory.

[0116] Based on this step, CytoEnergy identified four pathways for differentiation into terminal cells: Alpha, Beta, Delta, and Epsilon ( Figure 3 (C) in the figure, these four pathways represent the differentiation pathways of endocrine progenitor cells expressing low levels of Ngn3 into pancreatic endocrine cells that secrete different terminal hormones. It is worth mentioning that the Monocle (version 2.18.0) algorithm only identified one trajectory and could not correctly identify the four endocrine cell fate pathways ( Figure 3 (D) in the.

[0117] In addition, Gcg, Ins2, Sst, and Ghrl are known marker genes for Alpha, Beta, Delta, and Epsilon cells, respectively, and are specifically active in the corresponding cells ( Figure 4 (A1)-(A4) in Figure 2). Gcg is upregulated along the Alpha cell fate, and the local energy value of this gene also gradually increases along the Alpha path ( Figure 4 In addition, the activities of Ins2, Sst, and Ghrl gradually increase as they approach their associated terminal fate cells ( Figure 4 (B2)-(B4)), the expression levels of these genes were also upregulated ( Figure 4 The above results show that CytoEnergy technology not only captures the gene expression trend, but also takes into account the dynamic changes of the local gene network energy, providing a new perspective for in-depth characterization of the dynamic characteristics of core factors in cell fate transition.

[0118] Therefore, the present invention adopts the above-mentioned cell fate mapping research method based on CytoEnergy technology. Without relying on prior knowledge, this technology provides directional information for inferring cell fate transformation from the perspective of cell plasticity, and can accurately restore the dynamic change process of cell fate.

[0119] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A cell fate mapping research method based on CytoEnergy technology, characterized in that: The following steps are involved: Step S1, preprocessing single-cell transcriptome data; Step S2, constructing a local gene network based on the Hopfield neural network; Step S3, calculating the gene local network energy and cell energy, and normalizing the cell energy; Step S4, introducing a Gaussian mixture model to identify the potential state of cells by analyzing the difference in cell energy distribution, and clustering cells in the same differentiation potential state; Step S5, identifying cells in different terminal fate states; Step S6, inferring cell lineage trajectories; In step S5, cells in different terminal fate states are identified, and the specific operations are as follows: First, cell populations with different differentiation potentials were identified; Next, all cells were ranked according to their relative positions in the entire data set; When cells are arranged in chronological order from low to high in the actual biological process, the CytoEnergy technology defines the relative position ranking score of cells in the process of decreasing plasticity by Z'=1-E', where E' represents the normalized energy level of the cell; conversely, in the biological process of increasing cell plasticity, the position ranking of cells is determined by Z'=E'; CytoEnergy introduces a Gaussian distribution model, whose probability density function is defined as follows: Where μ represents Gaussian distribution The expected value of ; σ is the standard deviation, respectively expressed as: μ=median(Z′) σ=median(|Z′-μ|) Calculate Gaussian distribution The theta quantile of By default: θ = 85% where the position ranking score Z' exceeds the threshold cells, i.e. is considered to be in a potential terminal fate state; in biological processes with reduced plasticity, a low-potential state and Cells that are in a terminal fate state are considered to be in a terminal fate state; in contrast, in biological processes with increased plasticity, cells that simultaneously meet the high-potential state and Cells with a terminal fate are considered to be in a terminal fate state, and the pseudotime value of cells in a terminal fate state is set to 1.

2. A cell fate mapping research method based on CytoEnergy technology according to claim 1, characterized in that: In step S1, the single-cell transcriptome data is preprocessed, and the specific operations are as follows: Given a single-cell expression matrix containing N cells and M1 genes First, the single-cell transcriptome data was log-transformed using the NormalizeData function, and then standardized using the ScaleData function. Then, the highly variable genes were screened using the feature selection method. Finally, the highly variable genes were identified using the FindVariableFeatures function. Processed dataset contains N cells and M highly variable genes, and the state of cell n is represented by the vector express; Represents the state value of gene i in cell n. The specific calculation formula is as follows: in, Represents the normalized expression level of gene i in cell n.

3. The method for cell fate mapping based on CytoEnergy technology according to claim 2, characterized in that: In step S2, the nodes in the gene local network represent genes, and the interactions between genes in the network are specified by their adjacency matrix W. The matrix elements w ij Represents the relationship between genes i and j: in, represents the state value of gene j in cell n; p represents the number of stable states.

4. The method for cell fate mapping based on CytoEnergy technology according to claim 3, characterized in that: In step S3, the calculation formula of the gene local network energy is as follows: in, is the energy value of the local network centered on gene i in cell n; N(i) represents the neighborhood of node i in the network.

5. The method for cell fate mapping based on CytoEnergy technology according to claim 4, characterized in that: In step S3, the cell energy consists of the local network energy levels of all genes contained in the cell and is defined as: Among them, E n is the energy value of cell n; M represents the total number of highly variable genes.

6. The method for cell fate mapping based on CytoEnergy technology according to claim 5, characterized in that: In step S3, the normalization process is: in, and Respectively represent the minimum and maximum energy values ​​of all cells; E' n is the normalized cellular energy.

7. The method for cell fate mapping based on CytoEnergy technology according to claim 6, characterized in that: In step S4, a Gaussian mixture model is introduced to identify the potential state of cells by analyzing the energy distribution difference of cells, and cells in the same differentiation potential state are clustered. The specific operation is as follows: For a cell population S containing N cells, the normalized cell energy E' n∈N To quantify the differentiation potential status of each cell; The cells are sorted according to their energy and the E' of all cells is calculated using a Gaussian mixture model. n∈N The values ​​are fitted; Among them, the Gaussian mixture model is as follows: in, represents the kth Gaussian distribution with a given mean μ k and variance The probability density function under the condition λ k >0 is the mixing weight of the k-th Gaussian distribution, satisfying The number of components K of the Gaussian distribution is determined by the Bayesian Information Criterion BIC, which is defined as: Here, v is the number of estimated parameters.

8. The method for cell fate mapping based on CytoEnergy technology according to claim 7, characterized in that: In step S6, the cell lineage trajectory is inferred, and the specific operations are as follows: The cell fate of the remaining cells is measured by the combined similarity index CSI, which is defined as follows: CSI=(1-GD lt )*PC lt Among them, GD lt represents the Gower distance between the remaining cell l and the terminal state t, and is calculated as follows: in, represents the normalized expression level of gene i in cell l; is the average normalized expression level of gene i in all cells in terminal state t; d refers to the total number of genes in cell l; PC lt represents the Pearson correlation coefficient between the remaining cells l and the terminal state t, which is defined as follows: According to the CSI maximization principle, each remaining cell l is assigned to the differentiation path P represented by the corresponding terminal state t t As shown below: Where T represents the total number of terminal cell fate states; L refers to the total number of remaining cells, c l Represents the remaining cells l; c t represents the terminal state t; the pseudo-time value of cell l is determined by its position score Z', which reflects the relative position of the cell in the fate trajectory; The local weighted scatter point smoothing regression method was used, and the span parameter was set to 0.75 to fit the expression trend of genes along the lineage path; A span parameter of 0.95 was set to analyze the changing trend of the local network energy of genes along the lineage trajectory.

Images