Method for identifying cell fate transformation approach by using gene regulation network information

Abstract

The invention provides a method for identifying a cell fate transformation pathway by using gene regulation network information. The method comprises the following steps: constructing a co-expression gene network; calculating a cell fate probability according to the co-expression gene network; reconstructing a development track of cell fate transformation according to the cell fate probability; judging the activation state of the regulators according to the co-expression gene network; constructing a transcription factor feedback loop according to the co-expression gene network; and according to the development trajectory of the cell fate transformation, the activation state of the regulator and a transcription factor feedback loop, identifying a cell fate transformation path. According to the method, gene regulation network information for driving cell fate transformation can be captured, and a cell fate transformation path can be decoded according to the gene regulation network information.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a method for identifying a cell fate transition pathway using gene regulatory network information. Background Art

[0002] Cell fate transition refers to the dynamic process of cells changing from one function or differentiation state to another, including but not limited to cell differentiation, somatic cell reprogramming, dedifferentiation and transdifferentiation. Identifying and characterizing the complex pathways of cell fate transition and deciphering the intrinsic molecular mechanisms that drive these processes can help scientists understand and describe the process of cell fate transition more accurately, and provide a reference for identifying new therapeutic targets and developing new therapeutic strategies during the occurrence and evolution of diseases. At present, the algorithms for studying cell fate transition are mainly divided into: trajectory inference algorithms, algorithms for identifying cell fate tendencies, and algorithms based on RNA velocity. Among them, the trajectory inference algorithm is also called the pseudo-time algorithm, which is mainly built based on single-cell transcriptome data. Its core goal is to discover the one-dimensional potential representation of the cell state to calculate pseudo-time, and at the same time, it can identify the topological structure of cells in the embedded space to characterize complex cell development pathways. The algorithm for identifying cell fate tendencies focuses on quantitatively evaluating the uncertainty and probability of cell fate decisions, thereby identifying the fate preference of early cells. The algorithm based on RNA velocity calculates the transcription, splicing and degradation rates by detecting the relative abundance of unspliced ​​RNA and mature RNA, thereby estimating the dynamic change rate of gene expression, and further predicts the changes in the transcriptional state of cells in the future.

[0003] However, these existing algorithms have not fully explored and utilized the gene regulatory network information that plays a decisive role in cell state and function for modeling, and it is difficult to capture the gene regulatory network information that drives cell fate transformation. Therefore, it is necessary to design a method that uses gene regulatory network information to identify cell fate transformation pathways. Summary of the invention

[0004] The purpose of the present invention is to provide a method for identifying cell fate transition pathways using gene regulatory network information, so as to capture the gene regulatory network information that drives cell fate transition and decipher the cell fate transition pathway based on the gene regulatory network information.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A method for identifying cell fate transition pathways using gene regulatory network information comprises the following steps:

[0007] Constructing a co-expressed gene network;

[0008] Calculate cell fate probabilities based on co-expressed gene networks;

[0009] Reconstruct the developmental trajectory of cell fate transitions based on cell fate probabilities;

[0010] Determine the activation status of regulators based on the co-expressed gene network;

[0011] constructing transcription factor feedback loops based on co-expressed gene networks;

[0012] Identify cell fate transition pathways based on their developmental trajectories, activation states of regulators, and transcription factor feedback loops.

[0013] Optionally, the co-expression gene network is an expression matrix; the specific steps of constructing the co-expression gene network include:

[0014] Construct a partial correlation coefficient matrix based on the partial correlation coefficients of genes;

[0015] The false positive regulatory relationships in the partial correlation coefficient matrix were removed through the gene regulatory network to obtain the co-expressed gene network.

[0016] Optionally, the cell fate probability is calculated based on the co-expressed gene network, and the specific steps include:

[0017] The transmission cost matrix was obtained by coupling adjacent time points in the co-expressed gene network through Bayesian Ridge regression;

[0018] The transmission cost matrix is ​​optimized by the optimal transmission algorithm to obtain the optimal transmission matrix; the optimal transmission algorithm is the entropy regularization algorithm; the calculation formula of the entropy regularization algorithm is: subjecttoΣ j R ij =p i ,Σ i R ij =q j ; where i is an element in the gene source distribution, j is an element in the gene target distribution, p is the gene source distribution, q is the gene target distribution, R ij is the transmission amount from i to j, C ij For R ij The cost, p i and q j are the probabilities of elements in p and q respectively, ∈ is the regularization parameter;

[0019] The comprehensive transmission matrix is ​​obtained by coupling non-adjacent time points in the co-expressed gene network according to the Markov chain and the optimal transmission matrix;

[0020] Calculate cell fate probabilities based on the integrated transmission matrix.

[0021] Optionally, the development trajectory of cell fate transition is reconstructed according to the cell fate probability, and the specific steps are: analyzing the cell fate probability according to a preset fate probability threshold to obtain a Boolean mask, and reconstructing the development trajectory of cell fate transition according to the Boolean mask.

[0022] Optionally, the activation state of the regulator is determined based on the co-expressed gene network, and the specific steps include:

[0023] Identify transcription factors and target genes in gene regulatory networks and use transcription factors and target genes as regulator sets;

[0024] Calculate partial correlations between transcription factors and target genes;

[0025] The activation state of the regulator is determined based on the value of the partial correlation.

[0026] Optionally, judging the activation state of the regulator according to the co-expressed gene network further includes: calculating the activity score of the regulator; the calculation formula of the activity score of the regulator is: in, is the regulator activity fraction, is the expression value of the transcription factor in the cell, is the average expression value of the target gene that is positively regulated in the regulator in the cell, is the average expression value of the negatively regulated target gene in the regulator in the cell, w1, w2 and w3 are and The weight of .

[0027] Optionally, the activation state of the regulator is determined according to the value of the partial correlation, and the specific steps are:

[0028] When the positive correlation value is greater than 0.5, it is judged that the transcription factor and the target gene are in an activation relationship, that is, the regulator is in an activated state;

[0029] When the negative correlation value is less than -0.5, it is judged that the transcription factor and the target gene are in an inhibitory relationship, that is, the regulator is in an inhibitory state.

[0030] Optionally, a transcription factor feedback loop is constructed based on the co-expressed gene network, and the specific steps include:

[0031] Select transcription factors that act as downstream regulatory targets in gene regulatory networks;

[0032] Calculate the partial correlation between transcription factors and determine the functional relationship between transcription factors based on the partial correlation;

[0033] A directed graph is constructed based on the interaction relationship between transcription factors, and the transcription factor feedback loop is identified based on the directed graph.

[0034] Optionally, constructing a transcription factor feedback loop according to the co-expressed gene network further includes: calculating an activity score of the transcription factor feedback loop; the specific steps include:

[0035] Calculate the normalized mean of transcription factors in the transcription factor feedback loop; the calculation formula is: Among them, mean exp is the normalized mean of transcription factors, X i is the scaled expression value of the ith gene in the transcription factor feedback loop, and N is the number of genes in the transcription factor feedback loop;

[0036] The standard deviation was calculated based on the normalized mean of the transcription factors; the calculation formula was: Among them, std is the standard deviation;

[0037] The coefficient of variation was calculated based on the normalized mean and standard deviation of the transcription factors; the calculation formula is: Wherein, CV is the coefficient of variation;

[0038] The transcription factor feedback loop activity score was calculated based on the normalized mean and coefficient of variation of the transcription factor; the calculation formula is: TFFloop score =mean exp ×(1-CV); where TFFloop score is the transcription factor feedback loop activity fraction.

[0039] A system for identifying cell fate transition pathways using gene regulatory network information, comprising:

[0040] Network construction module: used to construct co-expression gene networks;

[0041] Data processing module: used to calculate cell fate probability based on the co-expressed gene network and reconstruct the development trajectory of cell fate transition based on the cell fate probability;

[0042] State judgment module: used to judge the activation state of regulators based on the co-expressed gene network;

[0043] Circuit identification module: used to construct transcription factor feedback loops based on co-expressed gene networks;

[0044] Cell recognition module: used to identify cell fate transition pathways based on the developmental trajectory of cell fate transition, the activation state of regulators, and the transcription factor feedback loop.

[0045] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects: the method provided by the present invention for identifying cell fate transition pathways using gene regulatory network information, the method comprising: constructing a co-expression gene network; calculating the cell fate probability based on the co-expression gene network; reconstructing the development trajectory of cell fate transition based on the cell fate probability; judging the activation state of the regulator based on the co-expression gene network; constructing a transcription factor feedback loop based on the co-expression gene network; identifying the cell fate transition pathway based on the development trajectory of the cell fate transition, the activation state of the regulator and the transcription factor feedback loop. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0047] Figure 1 This is a flow chart of a method for identifying cell fate transition pathways using gene regulatory network information according to an embodiment of the present invention;

[0048] Figure 2 This is a flow chart of calculating cell fate probability according to an embodiment of the present invention;

[0049] Figure 3 A schematic diagram of the probability of selecting cell fate according to an embodiment of the present invention;

[0050] Figure 4 A schematic diagram of a cell differentiation pathway selected according to an embodiment of the present invention;

[0051] Figure 5 A schematic diagram of the cell fate transformation development trajectory according to an embodiment of the present invention;

[0052] Figure 6 This is a flow chart of determining the activation state of a regulator according to an embodiment of the present invention;

[0053] Figure 7 This is a schematic diagram of constructing a transcription factor feedback loop according to an embodiment of the present invention. DETAILED DESCRIPTION

[0054] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0055] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and understandable, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0056] like Figure 1 As shown, an embodiment of the present invention provides a method for identifying a cell fate transition pathway using gene regulatory network information, comprising the following steps:

[0057] Step 100: Construct a co-expressed gene network.

[0058] Specifically, the co-expression gene network is an expression matrix based on single-cell transcriptome data. The rows of the matrix are gene names, the columns are cell labels, and the values ​​are the expression levels of genes in cells. The expression level is the logarithm of the gene reading value with base 2. The specific steps of constructing a co-expression gene network include:

[0059] First, calculate a variance-covariance matrix C X , X represents the set of all genes, C X The element c ij Represents x i and x j The covariance (x i and x j Represents the vector composed of genes and the expression level of genes in each cell respectively), and its calculation formula is: C X = Cov(x i ,x j ), where the diagonal elements c ii is x i The variance of Cov() is the operation to calculate the covariance matrix.

[0060] Next, calculate C X The inverse matrix of: inverse variance-covariance matrix D X ,Right now According to D X Calculate the partial correlation coefficient r ij|S , the calculation formula is: Among them, d ij d ii and d jj D X The elements in the i-th row and j-th column, the i-th row and i-th column, and the j-th row and j-th column in .

[0061] It should be noted that through the inverse variance-covariance matrix D X The elements in can directly calculate the correlation between two variables controlling other variables. This calculation method is relatively simple and can avoid tedious recursive calculations directly through matrix operations. It is especially suitable for high-dimensional data such as gene-cell expression matrix.

[0062] Again, by calculating the partial correlation coefficient between any two genes, we can get a The partial correlation coefficient matrix G X , the elements on the diagonal are null values, and this matrix is ​​also called the co-expression matrix. Where N is the number of genes and the number of matrix rows is The number of columns is 3, and these three columns represent gene i, gene j, and the partial correlation coefficient between gene i and gene j respectively.

[0063] Finally, the transcription factors in the regulatory DNA sequences (promoters and enhancers) of the open chromatin sites are scanned through the gene regulatory network (base gene regulatory network, base GRN), and the regulatory candidate genes are identified in combination with the motifs to remove the false positive regulatory relationships in the partial correlation coefficient matrix and obtain the co-expression gene network. The specific implementation method is: use a script customized in Python language to remove the regulatory relationships that do not exist in the partial correlation coefficient matrix according to the regulatory relationships of the base GRN, and adjust the data order of the three columns of the partial correlation coefficient matrix. The adjusted data order is: the first column is the transcription factor, the second round is the target gene, and the third column is still the partial correlation coefficient. This adjustment method narrows the scope of identifying regulatory candidate genes and helps to define the directionality of regulatory edges in the base GRN.

[0064] Step 200: Calculate cell fate probability based on the co-expressed gene network. The specific steps are as follows Figure 2 As shown, including:

[0065] Step 201: The transmission cost matrix is ​​obtained by coupling adjacent time points in the co-expression gene network through the Bayesian Ridge regression method. The specific implementation process is as follows:

[0066] First, the regulatory relationship between the target gene and its inferred upstream transcription factor was analyzed by Bayesian Ridge regression method, and the target gene was used as the dependent variable and the upstream transcription factor as the independent variable to obtain the Bayesian Ridge regression result, which is expressed as: Where y is the target gene, It indicates that the model assumes that the observed value follows a normal distribution, μ is the mean of the normal distribution and represents the expected linear prediction value, α is the intercept term, X is the upstream transcription factor of the target gene, β is the coefficient vector, and ∈ is the error term. Xβ can be expanded as: Among them, M is the number of upstream transcription factors inferred by the target gene, β m is the coefficient vector of the mth upstream transcription factor, x m is the mth upstream transcription factor. And the coefficient is used to further divide β into positive regulation (pos) transcription factors and negative regulation (neg) transcription factors, namely: k+l=M; where xpos is the expression of a positively regulated transcription factor, and its coefficient γ pos >0,x neg is the expression of a negatively regulated transcription factor, and its coefficient δ ng <0, pos and neg are variables, which represent the posth positively regulated transcription factor and the negth negatively regulated transcription factor, respectively. k is the number of positively regulated transcription factors, l is the number of negatively regulated transcription factors, and the sum of k and l is equal to the total number of upstream transcription factors M.

[0067] Then, two cell sets at adjacent time points t and t+1 are defined, and the unfavorable factors of the change in expression pattern of gene y between the two cell sets at t and t+1 are taken as the cost of the change in expression pattern of gene y between the two cell sets at t and t+1. These unfavorable factors include: the expression change of target gene Δy a,b , and Δy a,b The expression changes of upstream positive regulatory transcription factors with opposite signs Right now and Δy a,b The expression changes of upstream negative regulatory transcription factors with the same sign Right now 0 and.

[0068] Next, the transmission cost matrix of the expression pattern transition of gene y between the two cell sets at t and t+1 is calculated. The calculation formula is:

[0069]

[0070] Among them, a represents a cell at time point t, b represents a cell at time point t+1, represents the transmission cost of gene y from the expression pattern in cell a to the expression pattern in cell b, Δy a,b represents the difference in the expression level of gene y in cell a and cell b, represents the difference in expression level of the mth upstream transcription factor of gene y in cells a and b, β m is the regression coefficient corresponding to the mth upstream transcription factor of gene y, ads represents the absolute value, It is the condition for the calculation formula to be valid.

[0071] The transfer cost between any two cells is calculated by the cost of each gene between the two cells and the Euclidean distance of the cells in the two-dimensional principal component analysis space. The calculation formula is:

[0072]

[0073] Among them, cost a,bis the total cost of transforming cell a into cell b, sum(gene_cost a,b ) is the sum of the transmission costs of all genes, Represents the Euclidean distance between cell a and cell b in the two-dimensional principal component analysis space.

[0074] Finally, the transmission cost matrix of any two cells between the two time points t and t+1 is obtained based on the transfer cost between the two cells.

[0075] Step 202: Optimize the transmission cost matrix through the optimal transmission algorithm to obtain the optimal transmission matrix. The optimal transmission algorithm is the entropy regularization algorithm, and the calculation formula of the entropy regularization algorithm is:

[0076]

[0077] subjecttoΣ j R ij =p i ,Σ i R ij =q j ;

[0078] Where i is an element in the gene source distribution, j is an element in the gene target distribution, p is the gene source distribution, q is the gene target distribution, R ij is the transmission amount from i to j, C ij For R ij The cost, p i and q j are the probabilities of the elements in p and q respectively, and ∈ is the regularization parameter used to control the influence of the entropy regularization term. The entropy regularization algorithm can be used to obtain the optimal transmission matrix for any adjacent time points, and make the optimization problem smoother and easier to solve.

[0079] It should be noted that the known conditions in the optimal transmission algorithm are source distribution, target distribution and transmission cost. The transmission cost matrix is ​​the optimal transmission mode that minimizes the total transmission cost. The source distribution and target distribution are obtained through the growth rate. The growth rate is obtained by calculating the characteristics of the proliferation and apoptosis-related gene sets. The growth rate can intuitively indicate the speed of cell growth and is used to measure the potential of cells to produce offspring. The proliferation-related gene set is mainly involved in cell cycle-related processes, while the apoptosis-related gene set is mainly involved in cell autophagy-related processes. The calculation process of the growth rate is as follows:

[0080] First, the scores of the cell cycle and autophagy gene sets in each cell are calculated respectively; for any cell i, a gene set G is given, and the Z score of each gene j in the gene set is calculated using the following formula: Among them, X ijis the expression level of the jth gene in cell i, μ j and σ j are the mean expression and standard deviation of gene j, Z ij is the fraction of the jth gene in cell i.

[0081] The Z score is then truncated to [-5,5] to reduce the impact of extreme values. The calculation formula is: Z' ij =min(max(Z ij ,-5),5); where Z' ij is the Z score after truncation. Then according to Z' ij Calculate the average value of the gene set score corresponding to each cell using the following formula: Among them, S i is the final score of the cell cycle and autophagy gene set in the ith cell, and |G| is the number of genes in the gene set.

[0082] The final scores of the cell cycle and autophagy gene sets were then further converted to birth rate and death rate using the Logistic function to map the original gene set Z score to a biologically meaningful range. The calculation formula is: Among them, Rate represents the birth rate or death rate, and e is the base of the natural logarithm.

[0083] Finally, the cell growth rate (Growth rate) is calculated based on the birth rate and death rate. The calculation formula is: Growth rate = e Birthrate-Deathrate ; and through a certain time point t n and the next time point t n+1 The growth rates of cells are obtained by source distribution and target distribution.

[0084] Step 203: coupling non-adjacent time points in the co-expression gene network according to the Markov chain and the optimal transmission matrix to obtain a comprehensive transmission matrix.

[0085] Specifically, for any adjacent time node {t n ,t n+1}, its optimal transmission matrix is The next time node after the adjacent time node {t n+1 ,t n+2}, its optimal transmission matrix is Using the principle of Markov chain and Calculate and get non-adjacent time nodes {t n ,t n+2 The optimal transmission matrix of The calculation formula is: Then, by continuously right-multiplying the optimal transmission matrix of other adjacent time nodes, we can obtain N-1 With the last time node t N The optimal transmission matrix between them is obtained, and the comprehensive transmission matrix between all cells and the cells at the last time node is obtained through the optimal transmission matrix between any time node and the last time node.

[0086] Step 204: Calculate cell fate probability based on the integrated transmission matrix.

[0087] Specifically, the integrated transmission matrix is ​​normalized so that the sum of each row is 1, and the probability of each cell eventually transforming into the cell at the last time node is obtained. The cell at the last time node is taken as the final fate of the cell, and the columns of the integrated transmission matrix are summed according to the cell type to which the final fate of the cell belongs to obtain the cell fate probability. The schematic diagram of calculating the cell fate probability is shown in Figure 3 shown.

[0088] Step 300: Reconstruct the development trajectory of cell fate transformation according to cell fate probability. The specific steps are: analyzing the cell fate probability according to a preset fate probability threshold to obtain a Boolean mask, identifying the differentiation path of the cell according to the Boolean mask, and reconstructing the development trajectory of the cell fate transformation according to the differentiation path of the cell.

[0089] Specifically, the preset fate probability threshold is obtained by screening the fate probability distribution of all cells through a specific quantile, and the specific quantile in this embodiment is 99%. Figure 4 As shown in FIG. 1 , when the cell fate probability is higher than the preset fate probability threshold, the differentiation path of the cell is identified and stored. The dark part in the figure is the cell with a fate probability higher than the preset fate probability threshold. Figure 5 As shown, the differentiation pathways of multiple identified cells are combined to obtain the development trajectory of cell fate transformation (the dark part in the figure).

[0090] More specifically, the preset fate probability threshold adopts a dynamic adjustment strategy, divides the cell timeline into multiple sections, and then performs independent calculations in each section to obtain the quantile threshold of each section. The Boolean mask of each cell is obtained through the quantile threshold, and the quantile threshold is monotonically increased through the Boolean mask and the numpy.maximum.accumulate function of python. As long as the cell fate probability exceeds the corresponding dynamic threshold at any time point, the differentiation path of the cell is identified and stored.

[0091] It can be understood that the use of dynamic fate probability thresholds reflects the dynamic changes of cells during the differentiation process, improves the accuracy of identifying the differentiation path of cells, makes the development trajectory of reconstructed cell fate transformation more consistent with the real development trajectory, and greatly improves the accuracy of the reconstruction process.

[0092] Step 400: Determine the activation state of the regulator based on the co-expressed gene network. Figure 6 As shown, including:

[0093] Step 401: Identify transcription factors and target genes in a gene regulatory network, and use the transcription factors and target genes as a regulator set.

[0094] Specifically, the expression of genes in the regulatory subgroup is directly or indirectly regulated by the same transcription factor.

[0095] Step 402: Calculate the partial correlation between the transcription factor and the target gene.

[0096] Step 403: Determine the activation state of the regulator based on the value of the partial correlation.

[0097] Specifically, the specific steps for determining the activation state of the regulator are:

[0098] When the positive correlation value is greater than 0.5, it is judged that the transcription factor and the target gene are in an activation relationship, that is, the regulator is in an activated state; when the negative correlation value is less than -0.5, it is judged that the transcription factor and the target gene are in an inhibition relationship, that is, the regulator is in an inhibition state; the partial correlation value in other cases is adjusted to 0, indicating that there is no obvious regulatory relationship between the transcription factor and the target gene.

[0099] It should be noted that the target genes in the regulators are divided into two types: activated and inhibited by transcription factors through the judgment results of the activation state of the regulators. This classification method can identify regulators closely related to specific cell types. The specific implementation process is: each cell type is compared with all other cell types in turn, and the t-test is used to evaluate the expression difference of each regulator in different cells. When the p value in the t-test is less than 0.05, the regulator is considered to be significantly different from the cell. When a regulator is significantly different in comparison with all other cell types, it is considered to be a regulator closely related to a specific cell type.

[0100] More specifically, judging the activation state of the regulator based on the co-expressed gene network also includes: calculating the regulator activity score. The regulator activity score is used to evaluate the activity and functional importance of a specific transcription factor under a specific cell state or condition. The regulator activity score is obtained by weighting the contribution of transcription factor expression, activated target gene expression, and inhibited target gene expression, and the calculation formula is:

[0101]

[0102] in, is the activity fraction of regulator r in cell i, is the expression value of transcription factor in cell i, is the average expression value of the target gene positively regulated in regulator r in cell i, is the average expression value of the negatively regulated target gene in regulator r in cell i, w1, w2 and w3 are and In this embodiment, the above expression values ​​are scaled to between 0 and 1, w1=0.5, w2=w3=0.25.

[0103] Step 500: Construct a transcription factor feedback loop based on the co-expressed gene network. The specific steps are as follows: First, select transcription factors that can serve as downstream regulatory targets of other transcription factors to narrow the scope of analysis and form a set of transcription factors that may form a feedback loop. Next, determine whether there is an interaction between the transcription factors based on the partial correlation between the selected transcription factors. When the partial correlation is greater than 0.5, it is considered that there is a positive regulatory interaction between the transcription factors. Finally, each transcription factor is regarded as a node, and the positive regulatory interaction relationship between them is regarded as a directed edge to construct a directed graph, and the transcription factor feedback loop is constructed through a simple loop in the directed graph.

[0104] Specifically, constructing a transcription factor feedback loop according to the co-expressed gene network also includes: calculating the transcription factor feedback loop activity score; the specific steps include:

[0105] Calculate the normalized mean of all transcription factors in a transcription factor feedback loop in a cell; the calculation formula is: Among them, mean exp is the normalized mean of transcription factors, X i is the scaled expression value of the ith gene in the transcription factor feedback loop, and N is the number of genes in the transcription factor feedback loop. The normalized mean reflects the overall expression level of the transcription factor feedback loop in different cells.

[0106] The standard deviation was calculated based on the normalized mean and biased estimates of the transcription factors; the calculation formula is: Where std is the standard deviation, which is used to assess the consistency or variability of gene expression within a transcription factor feedback loop.

[0107] The coefficient of variation was calculated based on the normalized mean and standard deviation of the transcription factors; the calculation formula is: Among them, CV is the coefficient of variation.

[0108] The transcription factor feedback loop activity score was calculated based on the normalized mean and coefficient of variation of the transcription factor; the calculation formula is: TFFloop score =mean exp ×(1-CV); where TFFloop score is the activity score of the transcription factor feedback loop. This activity score is used to balance the overall expression level of the transcription factor feedback loop and the consistency of gene expression within the transcription factor feedback loop. The higher the activity score, the more active the transcription factor feedback loop is and the more consistent the gene expression is. Figure 7 As shown, the curve represents the change in the activity score of the transcription factor feedback loop, the gray circle represents the transcription factor, the dotted line with an arrow represents a positive regulatory interaction relationship that has not yet been established, and the solid line with an arrow represents an established positive regulatory interaction relationship.

[0109] Step 600: Identify cell fate transition pathways based on the developmental trajectory of cell fate transition, the activation state of regulators, and the transcription factor feedback loop.

[0110] A system for identifying cell fate transition pathways using gene regulatory network information, comprising:

[0111] Network construction module: used to construct co-expression gene networks;

[0112] Data processing module: used to calculate cell fate probability based on the co-expressed gene network and reconstruct the development trajectory of cell fate transition based on the cell fate probability;

[0113] State judgment module: used to judge the activation state of regulators based on the co-expressed gene network;

[0114] Circuit identification module: used to construct transcription factor feedback loops based on co-expressed gene networks;

[0115] Cell recognition module: used to identify cell fate transition pathways based on the developmental trajectory of cell fate transition, the activation state of regulators, and the transcription factor feedback loop.

[0116] The beneficial effects of the present invention are as follows:

[0117] 1) It achieves precise capture of regulators, greatly improving the depth and accuracy of identifying cell fate transition pathways;

[0118] 2) Combined with the dynamic changes of cells during differentiation, the accuracy of identifying the differentiation path of cells is improved;

[0119] 3) The capture of gene regulatory network information that drives cell fate transition was achieved, and the cell fate transition pathway was deciphered based on the gene regulatory network information.

[0120] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0121] The present invention uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only used to help understand the method and core ideas of the present invention. At the same time, for those skilled in the art, according to the ideas of the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for identifying cell fate transition pathways using gene regulatory network information, characterized in that: The steps include: Constructing a co-expressed gene network; calculating cell fate probabilities based on the co-expressed gene network; reconstructing the developmental trajectory of cell fate transitions according to the cell fate probabilities; Determining the activation state of the regulator according to the co-expressed gene network; constructing a transcription factor feedback loop based on the co-expressed gene network; The cell fate transition pathway is identified based on the developmental trajectory of the cell fate transition, the activation state of the regulator and the transcription factor feedback loop.

2. The method for identifying cell fate transition pathways using gene regulatory network information according to claim 1, characterized in that: The co-expression gene network is an expression matrix; the specific steps of constructing the co-expression gene network include: Construct a partial correlation coefficient matrix based on the partial correlation coefficients of genes; The co-expression gene network is obtained by removing the false positive regulatory relationships in the partial correlation coefficient matrix through the gene regulatory network.

3. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 1, characterized in that: Calculating cell fate probability according to the co-expressed gene network, the specific steps include: The transmission cost matrix is ​​obtained by coupling adjacent time points in the co-expression gene network through Bayesian Ridge regression method; The transmission cost matrix is ​​optimized by an optimal transmission algorithm to obtain an optimal transmission matrix; the optimal transmission algorithm is an entropy regularization algorithm; the calculation formula of the entropy regularization algorithm is: Where i is an element in the gene source distribution, j is an element in the gene target distribution, p is the gene source distribution, q is the gene target distribution, R ij is the transmission amount from i to j, C ij For R ij The cost, p i and q j are the probabilities of elements in p and q respectively, ∈ is the regularization parameter; coupling non-adjacent time points in the co-expression gene network according to the Markov chain and the optimal transmission matrix to obtain a comprehensive transmission matrix; The cell fate probabilities are calculated based on the integrated transmission matrix.

4. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 1, characterized in that: The development trajectory of cell fate transformation is reconstructed according to the cell fate probability, and the specific steps are: analyzing the cell fate probability according to a preset fate probability threshold to obtain a Boolean mask, and reconstructing the development trajectory of cell fate transformation according to the Boolean mask.

5. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 2, characterized in that: The activation state of the regulator is determined according to the co-expressed gene network, and the specific steps include: Identifying transcription factors and target genes in the gene regulatory network, and using the transcription factors and target genes as a regulator set; Calculating the partial correlation between the transcription factor and the target gene; The activation state of the regulator is determined according to the value of the partial correlation.

6. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 5, characterized in that: The method further comprises: determining the activation state of the regulator according to the co-expressed gene network, and calculating the activity score of the regulator; the calculation formula of the activity score of the regulator is: in, is the regulator activity fraction, is the expression value of the transcription factor in the cell, is the average expression value of the target gene that is positively regulated in the regulator in the cell, is the average expression value of the negatively regulated target gene in the regulator in the cell, w1, w2 and w3 are and The weight of .

7. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 5, characterized in that: The activation state of the regulator is judged according to the value of the partial correlation, and the specific steps are: When the positive correlation value is greater than 0.5, it is judged that the transcription factor and the target gene are in an activation relationship, that is, the regulator is in an activated state; When the negative correlation value is less than -0.5, it is judged that the transcription factor and the target gene are in an inhibitory relationship, that is, the regulator is in an inhibitory state.

8. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 2, characterized in that: Constructing a transcription factor feedback loop according to the co-expressed gene network, the specific steps include: Selecting a transcription factor that acts as a downstream regulatory target in the gene regulatory network; Calculating the partial correlations between the transcription factors, and determining the functional relationships between the transcription factors according to the partial correlations; A directed graph is constructed according to the transcription factors and the action relationships between the transcription factors, and the transcription factor feedback loop is identified according to the directed graph.

9. The method for identifying cell fate transition pathways using gene regulatory network information according to claim 8, characterized in that: Constructing a transcription factor feedback loop according to the co-expressed gene network also includes: calculating the activity score of the transcription factor feedback loop; the specific steps include: Calculate the normalized mean of the transcription factor in the transcription factor feedback loop; the calculation formula is: Among them, mean exp is the normalized mean of the transcription factor, X i is the scaled expression value of the ith gene in the transcription factor feedback loop, and N is the number of genes in the transcription factor feedback loop; The standard deviation was calculated based on the normalized mean of the transcription factor; the calculation formula was: Wherein, std is the standard deviation; The coefficient of variation is calculated based on the normalized mean value and the standard deviation of the transcription factor; the calculation formula is: Wherein, CV is the coefficient of variation; The transcription factor feedback loop activity score is calculated based on the normalized mean of the transcription factor and the coefficient of variation; the calculation formula is: TFFloop score =mean exp ×(1-CV); where TFFloop score is the activity score of the transcription factor feedback loop.

10. A system for identifying cell fate transition pathways using gene regulatory network information, applied to the method for identifying cell fate transition pathways using gene regulatory network information according to any one of claims 1 to 9, characterized in that: include: Network construction module: used to construct co-expression gene networks; Data processing module: used for calculating the cell fate probability according to the co-expressed gene network and reconstructing the development trajectory of cell fate transformation according to the cell fate probability; A state judgment module is used to judge the activation state of the regulator according to the co-expressed gene network; A loop identification module: used for constructing a transcription factor feedback loop according to the co-expressed gene network; Cell identification module: used to identify the cell fate transition pathway according to the development trajectory of the cell fate transition, the activation state of the regulator and the transcription factor feedback loop.

Images