A Method and System for Cell Pseudo-Time Trajectory Analysis Based on RNA Velocity and Gene Regulation Relationships

Abstract

A method and system for pseudo-time trajectory analysis of cells based on RNA velocity and gene regulatory relationships are disclosed. This method preprocesses gene expression to screen for potential regulatory relationships and establishes an RNA velocity model. The model parameters are iteratively optimized using a generalized EM algorithm. After integrating the optimized RNA velocity model across various genes, pseudo-time analysis and cell differentiation trajectory prediction are performed. This invention eliminates the dependence of existing technologies on spliced / unspliced ​​information. This significantly improves the model's fitting ability and extends the RNA velocity model to single-cell transcriptome datasets that cannot provide spliced / unspliced ​​information.

Description

Technical Field

[0001] This invention relates to a technology in the field of bioinformatics, specifically a method and system for analyzing cell pseudo-time trajectories based on the relationship between RNA velocity and gene regulation. Background Technology

[0002] RNA velocity describes the time derivative of gene expression. It predicts changes in gene expression and infers the fate of individual cells by modeling the relationship between the abundance of unspliced ​​(immature) and spliced ​​(mature) mRNAs. However, existing RNA velocity models can only be applied to scRNA-seq datasets with tags such as unspliced / spliced, and they cannot fit the dynamic curves of most genes well. Summary of the Invention

[0003] This invention addresses the shortcomings of existing technologies that rely solely on similarity to construct gene networks, neglecting gene regulatory relationships and failing to fit the temporal changes in gene expression. It proposes a cell pseudo-time trajectory analysis method and system based on RNA velocity and gene regulatory relationships. By modeling RNA velocity as a function of multiple transcription factor expression levels, it fully utilizes the transcriptional regulatory relationships between genes to establish dynamic equations for genes. Based on the expression levels of multiple transcription factors, it establishes dynamic change equations for target genes, thereby eliminating the dependence on splicing / unsplicing information in existing technologies. This significantly improves the model's fitting ability and extends the RNA velocity model to single-cell transcriptome datasets that cannot provide splicing / unsplicing information.

[0004] This invention is achieved through the following technical solution:

[0005] This invention relates to a cell pseudo-time trajectory analysis method based on RNA velocity and gene regulatory relationships. After preprocessing gene expression, potential regulatory relationships are screened out, and an RNA velocity model is established accordingly. The model parameters are iteratively optimized using the generalized EM algorithm. After integrating the RNA velocity optimization model into various genes, pseudo-time analysis and cell differentiation trajectory prediction are performed.

[0006] This invention relates to a cell pseudo-time trajectory analysis system based on RNA velocity and gene regulatory relationships to implement the above-mentioned method, comprising: a data preprocessing unit, a model optimization unit, and a downstream analysis unit, wherein: the data preprocessing unit preprocesses gene expression and screens for potential regulatory relationships to obtain the expression levels of each target gene and its corresponding transcription factor in each cell; the model optimization unit establishes an RNA velocity model for the expression level of each target gene and uses the EM algorithm to iteratively optimize the parameters in the model to obtain the optimal RNA velocity fitting model for the expression level of the target gene; the downstream analysis unit integrates the optimal RNA velocity fitting model for each target gene to predict cell development and differentiation trajectories. Attached Figure Description

[0007] Figure 1 This is a flowchart of the present invention;

[0008] Figure 2 This is a flowchart of an implementation example;

[0009] Figure 3 This is a schematic diagram illustrating the effect of the example;

[0010] In the figure: a) cell pseudo-time and cell trajectory predicted by TFvelo using the proposed method; b) cell pseudo-time and cell trajectory predicted by the existing scVelo; c) comparison of RNA velocity confidence; d) modeling effect on some genes; each row corresponds to the same gene, the first three columns show the dynamic model established in the state space by the proposed method, the distribution of gene expression on the UMAP plot, and the change of gene expression with pseudo-time; the last column shows the dynamic model modeled by the existing technology, and the absence of a fitting curve means that the existing technology has failed to establish any model. Detailed Implementation

[0011] like Figure 1 As shown in the figure, this embodiment relates to a method for predicting cell development and differentiation trajectories based on RNA velocity modeling, which includes the following steps:

[0012] Step 1) Preprocess the single-cell transcriptome data to obtain the expression levels of each target gene and its corresponding transcription factor in each cell, specifically including:

[0013] Step 1.1) Preprocess gene expression levels, specifically as follows:

[0014] Step 1.1.1) Filtering for low-expression genes. First, remove genes that are expressed in less than 2% of cells, and then retain only the top 2000 hypervariable genes, i.e., the 2000 genes with the largest expression variance.

[0015] Step 1.1.2) In each cell, the gene expression level is normalized, that is, the expression count of each gene is divided by the total gene expression count in each cell to obtain the expression level of each gene.

[0016] Step 1.1.3) Perform principal component analysis on the cells in the gene expression space, and calculate the 30 nearest neighbors of each cell based on the first 30 principal components.

[0017] Step 1.1.4) Smooth the gene expression data using the nearest neighbor algorithm. That is, for each cell, smooth the gene expression levels of each gene in the cell using the gene expression levels of its 30 neighbors to obtain preprocessed gene expression data.

[0018] Step 1.2) Extract transcription factor-target gene combinations with potential regulatory relationships. Specifically, download all labeled gene-gene transcriptional regulatory relationships from the ENCODE transcription factor-target gene dataset and the ChEA transcription factor-target gene dataset. If a set of transcription factor-target genes is labeled in at least one of the two datasets, add the transcription factor to the potential transcription factor array of the target gene to obtain the transcription factor expression data for each target gene.

[0019] Step 2) Establish an RNA velocity model based on gene regulatory relationships, and use the transcription factor expression data for each target gene obtained in Step 1 to optimize the model parameters using the generalized EM algorithm, specifically including:

[0020] Step 2.1) For each gene, establish a dynamic model of RNA velocity, specifically as follows:

[0021] Step 2.1.1) RNA velocity is modeled as a linear combination of transcription factors and the expression level of the gene, i.e. Wherein: the expression level of gene g is y g The expression levels of all transcription factors that have potential regulatory relationships with it are X. g , This represents the RNA velocity of gene g, W g W represents the weight of each transcription factor. g X g (t) represents the sum of the regulatory effects of transcription factors on gene g, γ g Indicates gene γ g The attenuation rate. W g and γ g The solution will be optimized in subsequent steps.

[0022] Step 2.1.2) The change in gene expression over time is modeled as a sinusoidal function, which satisfies the properties of high-order differentiability and a single peak, i.e., y g(t)=α g sin(2πt+θ g )+β g , t∈[0,1). where α g ,β g θ g These are all parameters in the function that characterizes gene expression over time, and will be optimized and solved in subsequent steps.

[0023] Step 2.2) W g X g Same as y g The joint space serves as the state space.

[0024] Step 2.3) In the state space, construct a loss function based on the distance from each cell to the model curve. Among them, the observed expression levels of transcription factors and target genes in cell c were: and

[0025] Step 2.4) The parameters are divided into cell time t and transcription factor weights W. g and model parameter set [α g ,β g θ g γ g The generalized EM algorithm is used for alternating optimization, with the optimization objective being to minimize the distance from each cell to the model curve. Specifically, this includes:

[0026] Step 2.4.1) Use grid search to assign time t to each cell, that is, locate the point on the model curve that is the closest to each cell, and regard this point as the target point of its corresponding cell. Steps 2.4.2) and 2.4.3) are dedicated to reducing the distance between each cell and its target point on the model curve in the state space.

[0027] Step 2.4.2) Optimize the weights W of the transcription factors using a linear regression method with inequality constraints, namely the Trustregion reflective algorithm. g This minimizes the loss function.

[0028] Step 2.4.3) Use the quasi-Newton method L-BFGS-B to test α g ,β g θ g γ g The value of is optimized to minimize the loss function.

[0029] Step 2.5) Optimizes the model parameters by iterating through steps 2.4.1) to 2.4.3) for 20 rounds. Since EM-type algorithms can only converge to local optima, this invention selects different initial parameter values ​​and applies the above optimization process accordingly. The result with the smallest error is retained, thus obtaining a fitting model for the dynamic changes of genes.

[0030] Step 3) Integrate the fitted models of dynamic changes for each gene obtained in Step 2, and perform pseudo-time analysis, cell trajectory prediction, and analysis of key transcriptional regulatory relationships, specifically including:

[0031] Step 3.1) The pseudo-time of cells is estimated using an improved method based on the literature (Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nature biotechnology. 2020; 38(12):1408-14.), specifically including:

[0032] Step 3.1.1) Calculate the cosine similarity between velocity and potential cell state transitions to obtain the transition matrix Π, where each element π in the matrix... ij The probability of cell i transitioning to cell j: Where: v i The RNA velocity of each gene in cell i, ≤ ij =s j -s i Measure the differences in gene expression between cell i and cell j.

[0033] Step 3.1.2) Perform matrix decomposition on the transformation matrix Π to obtain its eigenvalues ​​and eigenvectors. Identify the cells corresponding to the eigenvectors with an eigenvalue of 1 as the endpoint cells. Similarly, calculate the transpose Π of the transformation matrix. T The eigenvalues ​​and eigenvectors are used to identify the eigenvectors with a value of 1 as the starting cells.

[0034] Step 3.1.3) For each cell, calculate the average number of steps required to reach that cell from the root cell as the cell pseudo-time.

[0035] Step 3.2) Based on the predicted cell pseudo-time, predict the cell trajectory, specifically including:

[0036] Step 3.2.1) Based on the neighbor relationships between cells, use UMAP to reduce and visualize the data, thereby obtaining the pseudo-temporal distribution of cells on a two-dimensional UMAP map.

[0037] Step 3.2.2) Divide the UMAP map into 50x50 small grids. For each cell in the grid, calculate its average pseudo-time as the pseudo-time of that grid.

[0038] Step 3.2.3) Calculate the pseudo-time gradient for each grid on the UMAP plot and draw arrows based on these gradients. The arrows indicate the direction of cell differentiation.

[0039] Through specific experiments on single-cell transcriptome data, the effectiveness of this invention and its advantages over existing technologies were verified. Specifically, a dataset of erythrocyte developmental lines from the gastrula cytogenetic lineage of mouse embryos was used, containing the expression of 53,801 genes from 9,815 cells. This dataset has been widely used in previous RNA velocity models. After installing the scvelo package in Python, this data can be automatically downloaded and read using scvelo.datasets.gastrulation_erythroid().

[0040] like Figure 3 As shown, compared with the existing technology (scVelo), this method (TFvelo) can better reconstruct the pseudo-temporal distribution of cells and fit the cell differentiation trajectory (in reality, cells follow...). Figure 3 (Differentiation occurs in the order from the top left to the bottom right in cells a and 3b). Furthermore, velocity confidence reflects the consistency of RNA velocity among cells in similar states. Results are as follows... Figure 3 As shown in c, the RNA speed of this method has higher consistency.

[0041] The modeling effect of each gene is as follows: Figure 3 As shown in d. Because the splicing information used in existing techniques is highly noisy, it often fails to provide accurate modeling. In contrast, this method can fully utilize the transcriptional regulatory relationships between genes to better fit the actual state of the cell.

[0042] Compared with existing technologies, the model used in this invention can better fit the actual collected gene expression data, and its inferred cell differentiation trajectory and time situation are more consistent. The obtained RNA velocity distribution has higher neighbor consistency. In addition, since the RNA velocity modeling method proposed in this invention is free from dependence on splicing information in the data, the RNA velocity model can be extended to a wider range of single-cell transcriptome datasets, such as data measured by fluorescence in situ hybridization and other techniques.

[0043] The above-described specific implementations can be partially adjusted by those skilled in the art in different ways without departing from the principles and purpose of the present invention. The scope of protection of the present invention is defined by the claims and is not limited to the above-described specific implementations. All implementation schemes within the scope of the claims are bound by the present invention.

Claims

1. A method for predicting cell development and differentiation trajectories based on RNA velocity modeling, characterized in that, After preprocessing gene expression, potential regulatory relationships were screened, and an RNA velocity model was established accordingly. The model parameters were iteratively optimized using the generalized EM algorithm. After integrating the optimized RNA velocity model across various genes, pseudo-time analysis and cell differentiation trajectory prediction were performed, specifically including: Step 1) Preprocess the single-cell transcriptome data to obtain the expression levels of each target gene and its corresponding transcription factor in each cell, specifically including: Step 1.1) Preprocess gene expression levels; Step 1.2) Extract transcription factor-target gene combinations with potential regulatory relationships. Specifically, download all labeled gene transcriptional regulatory relationships from the ENCODE transcription factor-target gene dataset and the ChEA transcription factor-target gene dataset. If a set of transcription factor-target genes is labeled in at least one of the above two datasets, add the transcription factor to the potential transcription factor array of the target gene to obtain the transcription factor expression data for each target gene. Step 2) Establish an RNA velocity model based on gene regulatory relationships, and use the transcription factor expression data for each target gene obtained in Step 1) to optimize the model parameters using the generalized EM algorithm, specifically including: Step 2.1) For each gene, establish a dynamic model of RNA velocity, specifically including: Step 2.1.1) RNA velocity modeling is a linear combination of the expression levels of each transcription factor and the gene, i.e. Among them: genes The expression level is The expression levels of all transcription factors that have potential regulatory relationships with it are , That is, genes RNA speed, The weights of each transcription factor, For each transcription factor to gene The sum of regulatory effects, For genes The attenuation rate; and The solution will be optimized by subsequent steps; Step 2.1.2) Model the change in gene expression over time as a sinusoidal function, which satisfies the properties of high-order differentiability and a single peak, i.e. ;in: These are all parameters in the function that characterizes gene expression over time; Step 2.2) Set the parameters of the dynamic model , same The joint space serves as the state space; Step 2.3) In the state space, construct a loss function based on the distance from each cell to the model curve. Among them: the observed transcription factors and target genes in cells The expression levels in are respectively and ; Step 2.4) Divide the parameters into cell time. Weights of transcription factors and model parameter set [ The generalized EM algorithm is used for alternating optimization, with the optimization objective being to minimize the distance from each cell to the model curve; Step 2.5) By selecting different initial values ​​of parameters, iterate through Step 2.4) and select the result with the smallest error to obtain a fitting model for the dynamic changes of genes. Step 3) Integrate the fitted models of dynamic changes for each gene obtained in Step 2), and perform pseudo-time analysis, cell trajectory prediction, and analysis of key transcriptional regulatory relationships, specifically including: Step 3.1) Use an improved method to estimate the pseudo-time of cells; Step 3.2) Predict cell trajectories based on the predicted cell pseudo-time.

2. The method for predicting cell development and differentiation trajectories based on RNA velocity modeling according to claim 1, characterized in that, Step 1.1) specifically includes: Step 1.1.1) Filter low-expression genes; first remove genes that are expressed only in less than 2% of cells, and then retain only the top 2000 hypervariable genes, that is, the 2000 genes with the largest expression variance. Step 1.1.2) In each cell, the gene expression level is normalized, that is, the expression count of each gene is divided by the total gene expression count in each cell to obtain the expression level of each gene. Step 1.1.3) Perform principal component analysis on the cells in the gene expression space, and calculate the 30 nearest neighbors of each cell based on the first 30 principal components; Step 1.1.4) Use the nearest neighbor algorithm to smooth the gene expression data; that is, for each cell, use the gene expression levels of its 30 neighbors to smooth the gene expression levels of each gene in the cell to obtain preprocessed gene expression data.

3. The method for predicting cell development and differentiation trajectories based on RNA velocity modeling according to claim 1, characterized in that, Step 2.4) specifically includes: Step 2.4.1) Assign time t to each cell using a grid search method, that is, locate the point with the smallest distance from the model curve to each cell, and regard this point as the target point of its corresponding cell; Step 2.4.2) Optimize the weights of transcription factors using a linear regression method with inequality constraints, specifically the trust region reflection algorithm. This minimizes the loss function; Step 2.4.3) Use the quasi-Newton method to... The value of is optimized to minimize the loss function, thereby reducing the distance between each cell in the state space and its target point on the model curve.

4. The method for predicting cell development and differentiation trajectories based on RNA velocity modeling according to claim 1, characterized in that, Step 3.1) specifically includes: Step 3.1.1) Calculate the cosine similarity between velocity and potential cell state transitions to obtain the transition matrix. Each element in the matrix The probability of cell i transitioning to cell j: ,in: The RNA velocity of each gene within cell i. Measure the differences in gene expression between cell i and cell j; Step 3.1.2) Transformation matrix Perform matrix decomposition to obtain its eigenvalues ​​and eigenvectors; identify the cells corresponding to the eigenvectors with an eigenvalue of 1 as the endpoint cells; similarly, calculate the transpose of the transformation matrix. The eigenvalues ​​and eigenvectors are used to identify the eigenvectors whose eigenvalues ​​are 1 as the originating cells. Step 3.1.3) For each cell, calculate the average number of steps required to reach that cell from the root cell as the cell pseudo-time.

5. The method for predicting cell development and differentiation trajectories based on RNA velocity modeling according to claim 1, characterized in that, Step 3.2) specifically includes: Step 3.2.1) Based on the neighbor relationships between cells, use UMAP to reduce the dimensionality of the data and visualize it, thereby obtaining the pseudo-temporal distribution of cells on the two-dimensional UMAP map; Step 3.2.2) Divide the UMAP map into 50x50 small grids. For each cell in the grid, calculate the average pseudo time as the pseudo time of that grid. Step 3.2.3) Calculate the pseudo-time gradient for each grid on the UMAP plot and draw arrows based on these gradients. The arrows indicate the direction of cell differentiation.

6. A system for implementing the cell development and differentiation trajectory prediction method based on RNA velocity modeling as described in any one of claims 1-5, characterized in that, include: The system comprises a data preprocessing unit, a model optimization unit, and a downstream analysis unit. The data preprocessing unit preprocesses gene expression and screens for potential regulatory relationships to obtain the expression levels of each target gene and its corresponding transcription factor in each cell. The model optimization unit establishes RNA velocity models for the expression levels of each target gene and uses the EM algorithm to iteratively optimize the parameters in the models to obtain the optimal RNA velocity fitting model for the target gene expression levels. The downstream analysis unit integrates the optimal RNA velocity fitting models for each target gene to predict cell development and differentiation trajectories.

Images