A method for identifying cell fate transition pathways using gene regulatory network information

Abstract

The present invention provides a method for identifying cell fate transition pathways using gene regulatory network information. The method comprises: 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 based on the cell fate probabilities; determining the activation state of regulators based on the co-expressed gene network; constructing a transcription factor feedback loop based on the co-expressed gene network; and identifying cell fate transition pathways based on the developmental trajectory of cell fate transitions, the activation state of regulators, and the transcription factor feedback loop. This method can capture the gene regulatory network information that drives cell fate transitions and decipher the cell fate transition pathways based on the gene regulatory network information.

Description

Technical Field

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

[0002] Cell fate transition refers to the dynamic process by which cells change from one functional 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 transitions and deciphering the underlying molecular mechanisms driving these processes can help scientists more precisely understand and describe the processes of cell fate transition, providing insights for identifying new therapeutic targets and developing novel treatment strategies during disease development and progression. Currently, algorithms for studying cell fate transitions are primarily categorized into trajectory inference algorithms, algorithms for identifying cell fate propensities, and algorithms based on RNA velocity. Trajectory inference algorithms, also known as pseudotime algorithms, are primarily built based on single-cell transcriptome data. Their core goal is to discover a one-dimensional latent representation of cell states to calculate pseudotime. They can also identify the topological structure of cells in embedded space to characterize complex cellular developmental pathways. Algorithms for identifying cell fate propensities focus on quantitatively assessing the uncertainty and probability of cell fate decisions, thereby identifying fate preferences of early-stage cells. RNA velocity-based algorithms calculate transcription, splicing, and degradation rates by measuring the relative abundance of unspliced ​​and mature RNA, thereby estimating the dynamic rate of gene expression and further predicting changes in the transcriptional state of cells over time.

[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 have difficulty capturing the gene regulatory network information that drives cell fate transitions. Therefore, it is very necessary to design a method that uses gene regulatory network information to identify cell fate transition 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] Construct 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, regulator activation states, 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 co-expression gene network was obtained by removing the false positive regulatory relationships in the partial correlation coefficient matrix through the gene regulatory network.

[0016] Optionally, cell fate probabilities are calculated based on a co-expressed gene network, specifically comprising:

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

[0018] The optimal transmission matrix is ​​obtained by optimizing the transmission cost matrix through the optimal transmission algorithm; 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 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, specifically comprising the following steps: 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 them as regulator sets;

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

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

[0026] Optionally, determining the activation state of the regulator based on the co-expressed gene network further includes: calculating the regulator activity score; the calculation formula for the regulator activity score 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 based on the partial correlation value, and the specific steps are as follows:

[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 interaction 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 based on 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 the transcription factors 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 i-th 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: Where 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-expression 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 cell fate probabilities based on the co-expression gene network; reconstructing the development trajectory of cell fate transition based on the cell fate probabilities; judging the activation state of regulators based on the co-expression gene network; constructing a transcription factor feedback loop based on the co-expression gene network; and identifying cell fate transition pathways based on the development trajectory of cell fate transition, the activation state of regulators, 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 following briefly introduces the drawings required for use in the embodiments. 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 any creative work.

[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 Schematic diagram of the cell fate transition 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 Schematic diagram of the transcription factor feedback loop construction according to an embodiment of the present invention. DETAILED DESCRIPTION

[0054] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. 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 making creative efforts 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 easy to understand, 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 cell fate transition pathways 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 based on an expression matrix of 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 the genes in the cells. The expression levels are the logarithm of the gene read counts to the base 2. The specific steps to construct the 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 consisting of genes and gene expression levels in each cell respectively), and its calculation formula is: C X = Cov(x i ,x j ); among them, the elements c on the diagonal 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 at row i, column j, row i, column i, and row j, column j in .

[0061] It should be noted that 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 by directly performing matrix operations. It is particularly suitable for high-dimensional data such as gene-cell expression matrices.

[0062] Again, after 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.

[0063] Finally, the base gene regulatory network (base GRN) is used to scan the transcription factors within the regulatory DNA sequences (promoters and enhancers) of open chromatin sites, and the motifs are combined to identify regulatory candidate genes to remove false positive regulatory relationships in the partial correlation coefficient matrix and obtain a co-expression gene network. The specific implementation method is: using 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 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. Figure 2 As shown, including:

[0065] Step 201: Using Bayesian Ridge regression to couple adjacent time points in the co-expressed gene network, a transmission cost matrix is ​​obtained. 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. The target gene was used as the dependent variable and the upstream transcription factor as the independent variable. The Bayesian Ridge regression results were obtained and expressed as: Where y is the target gene, Indicates that the model assumes that the observed values ​​follow 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. 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 M of upstream transcription factors.

[0067] Then, define two cell sets at adjacent time points t and t+1, and take the unfavorable factors of the expression pattern transition of gene y between the two cell sets at t and t+1 as the cost of the transition of 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 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 as follows:

[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 its expression pattern in cell a to its expression pattern in cell b, Δy a,b represents the difference in expression levels of gene y in cell a and cell b, represents the difference in expression levels 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 using an optimal transmission algorithm to obtain an optimal transmission matrix. The optimal transmission algorithm is an 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 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. Using the entropy regularization algorithm, the optimal transmission matrix for any adjacent time points can be obtained, making the optimization problem smoother and easier to solve.

[0079] It should be noted that the known conditions in the optimal transmission algorithm are the source distribution, target distribution, and transmission cost. The transmission cost matrix is ​​the optimal transmission method 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 more 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 growth rate calculation process is as follows:

[0080] First, the scores of the cell cycle and autophagy gene sets in each cell are calculated separately; 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 level 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 i-th 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: Here, 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 from the source distribution and the 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 And the 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 multiplying the optimal transmission matrix of other adjacent time nodes on the right, we can get 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 probabilities 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 a cell at the last time node is obtained. The cell at the last time node is regarded 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 cell's final fate belongs to obtain the cell fate probability. The schematic diagram of calculating cell fate probability is shown in the figure below. Figure 3 shown.

[0088] Step 300: Reconstructing the developmental trajectory of cell fate transitions based on cell fate probabilities. Specifically, the steps include analyzing the cell fate probabilities based on a preset fate probability threshold to obtain a Boolean mask, identifying the cell differentiation pathway based on the Boolean mask, and reconstructing the developmental trajectory of cell fate transitions based on the cell differentiation pathway.

[0089] Specifically, the preset fate probability threshold is obtained by screening the fate probability distribution of all cells through a specific quantile. In this embodiment, the specific quantile is 99%. Figure 4 As shown in the figure, 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 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 thresholds are dynamically adjusted, dividing the cell's timeline into multiple segments. Within each segment, independent calculations are performed to determine the quantile threshold for each segment. This quantile threshold is used to generate a Boolean mask for each cell. This mask is then used with Python's numpy.maximum.accumulate function to monotonically increase the quantile threshold. Whenever a cell's fate probability exceeds the corresponding dynamic threshold at any point in time, the cell's differentiation path 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 transition 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 activation 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; in other cases, the partial correlation value 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 results of judging the activation status of regulators divide the target genes in the regulators into two types: activated by transcription factors and inhibited by transcription factors. This classification method can identify regulators closely associated with specific cell types. The specific implementation process is: each cell type is compared with all other cell types in turn, and a t-test is used to evaluate the expression differences 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 that cell type. When a regulator is significantly different in comparison with all other cell types, it is considered to be a regulator closely associated with a specific cell type.

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

[0101]

[0102] in, is the activity fraction of regulator r in cell i, is the expression value of the transcription factor in cell i, is the average expression value of the target gene positively regulated by 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 analysis scope 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. The transcription factor feedback loop is constructed through simple loops in the directed graph.

[0104] Specifically, constructing a transcription factor feedback loop based on 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 within a cell; 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 i-th gene in the transcription factor feedback loop, and N is the number of genes in the transcription factor feedback loop. This 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: Where 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-expression 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) Achieved precise capture of regulatory elements, greatly improving the depth and accuracy of identifying cell fate transition pathways;

[0118] 2) Incorporating the dynamic changes of cells during differentiation, the accuracy of identifying cell differentiation pathways is improved;

[0119] 3) The capture of gene regulatory network information that drives cell fate transitions was achieved, and the cell fate transition pathways were 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 intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed 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-expression gene network; the co-expression gene network is an expression matrix; the specific steps of constructing the co-expression gene network include: constructing a partial correlation coefficient matrix based on the partial correlation coefficients of genes; removing false positive regulatory relationships in the partial correlation coefficient matrix through a gene regulatory network to obtain the co-expression gene network; Calculating cell fate probability based on 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 using the Bayesian Ridge regression method. The calculation formula of the transmission cost matrix is: ;in, Indicates a time point A cell, Indicates a time point A cell, Indicates gene From cells The expression pattern in cells changes The transmission cost of the expression pattern in Indicates gene In cells and cells The difference in expression levels, Indicates gene No. Upstream transcription factors in cells and cells The difference in expression levels, It's genes No. The regression coefficient corresponding to the upstream transcription factor, Indicates taking the absolute value; 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, for arrive The amount of transmission, for The cost, and They are and The probability of an element in is the regularization parameter; coupling non-adjacent time points in the co-expressed gene network according to the Markov chain and the optimal transmission matrix to obtain a comprehensive transmission matrix; calculating the cell fate probability according to the integrated transmission matrix; Reconstructing the development trajectory of cell fate transition according to the cell fate probability; 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; Determining the activation state of the regulator based on the co-expressed gene network; 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 judged according to the value of the partial correlation; 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 activation 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; Constructing a transcription factor feedback loop based on 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 based on the partial correlations; constructing a directed graph according to the transcription factors and the interaction relationships between the transcription factors, and identifying the transcription factor feedback loop according to the directed graph; 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: Determining the activation state of the regulator based on the co-expressed gene network also includes: calculating the regulator activity score; the calculation formula of the regulator activity score 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, 、 They are 、 and The weight of .

3. The method of identifying cell fate transition pathways using gene regulatory network information according to claim 1, characterized in that: Constructing a transcription factor feedback loop based on the co-expressed gene network also includes: calculating the transcription factor feedback loop activity score; the specific steps include: Calculate the normalized mean of the transcription factor in the transcription factor feedback loop; the calculation formula is: ;in, is the normalized mean of the transcription factor, is the scaled expression value of the i-th 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 was calculated based on the normalized mean and standard deviation of the transcription factor; the calculation formula was: ; Wherein, CV is the coefficient of variation; The transcription factor feedback loop activity score is calculated based on the normalized mean value and the coefficient of variation of the transcription factor; the calculation formula is: ;in, is the activity score of the transcription factor feedback loop.

4. 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 3, characterized in that: include: Network construction module: used to construct co-expression gene networks; 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; State judgment module: used to judge the activation state of the regulator based on the co-expressed gene network; A circuit identification module is used to construct a transcription factor feedback circuit based on the co-expressed gene network; Cell identification module: used to identify 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.

Images