Dynamic gene regulatory network inference method based on optimal transmission and multi-omics prior fusion

Abstract

This invention discloses a dynamic gene regulatory network inference method based on optimal transport and multi-omics prior fusion. This method can accurately reconstruct the continuous evolution process of gene regulatory networks from discrete single-cell multi-omics data, effectively avoiding directional bias at key differentiation nodes. By fusing transcriptomic and chromatin accessibility information, the reliability of regulatory priors is significantly improved, enabling accurate differentiation of complex regulatory scenarios. The end-to-end graph neural network architecture allows regulatory inference and phenotypic prediction to mutually reinforce each other, and error signals can be backpropagated to optimize regulatory relationships, thereby obtaining a high-resolution, high-precision dynamic regulatory weight matrix. This method remains robust under sparse data conditions, adaptively balancing priors and exploration at key bifurcation states, significantly reducing false negative and false positive biases, and providing reliable computational support for analyzing dynamic processes such as development and disease.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and computational biology, specifically to a method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion. Background Technology

[0002] In the formation, regeneration, and disease progression of complex tissues, the dynamic synergistic effects of gene regulatory networks (GRNs) play a central role. Recent breakthroughs in single-cell sequencing technology have made it possible to simultaneously analyze the transcriptome and chromatin accessibility in the native cell state, providing a crucial data foundation for inferring the dynamic evolution of GRNs. However, reconstructing the true continuous dynamics of GRNs from discrete, sparse single-cell multi-omics snapshots still faces significant computational challenges. Existing methods suffer from fundamental flaws in handling the contradiction between temporal continuity and discrete sampling: actual collected data often only represent a few discrete time points, and techniques such as moving window smoothing and variational inference trajectory estimation essentially fail to address the problem of modeling the continuous evolution of cell states between adjacent time points, leading to a systematic amplification of directional deviations in regulatory relationships at key differentiation nodes. In the field of multi-omics fusion, current mainstream methods mostly integrate transcriptome and ATAC data by feature-level splicing or probability product, ignoring the biophysical hierarchy of chromatin openness as a regulatory "permission" signal. That is, the accessibility of transcription factor binding motifs should be regarded as a gating switch of their regulatory ability rather than a simple feature superposition. Therefore, it is impossible to accurately distinguish between scenarios such as "high expression of transcription factors but closed binding sites" and "low expression but open binding sites", which seriously limits the reliability of regulatory priors.

[0003] Furthermore, existing GRN inference methods mostly output coarse-grained trend curves of regulatory strength, lacking the ability to output specific regulatory matrices at each interpretable real time point. Simultaneously, GRN inference and phenotypic prediction remain in a separate two-stage paradigm, preventing the backpropagation of phenotypic prediction error signals to optimize the inference process, resulting in weak correlation between regulatory relationships and downstream phenotypic tasks. Finally, existing models either rely entirely on data-driven approaches, making them prone to overfitting in data-sparse situations, or they over-rely on fixed priors, lacking the ability to explore unannotated regulatory relationships. They cannot adaptively balance the constraints of prior knowledge with data-driven exploration when cells transition from a steady state to a critical bifurcation state, thus producing false negatives or false positives in the prediction of key regulatory events. In summary, there is an urgent need for a GRN inference method that can integrate multi-omics data, reconstruct continuous evolutionary processes from discrete time points, output high-resolution regulatory weight matrices, and simultaneously predict cell phenotypes. Summary of the Invention

[0004] To address the technical problems mentioned above, this invention aims to provide a time-resolved gene regulatory network inference method based on optimal transport interpolation, RNA rate-driven ordinary differential equations, multi-omics prior fusion, and graph neural networks. This method enables the reconstruction of continuous gene expression dynamics from discrete time-point single-cell data; inference of the time-varying regulatory weights between transcription factors and target genes; integration of chromatin openness prior constraints provided by ATAC-seq; and simultaneous output of cell phenotype predictions, thereby improving the model's interpretability and application value.

[0005] To achieve the above objectives, this invention provides a method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion, comprising the following steps: S1. Based on single-cell transcriptome data and RNA rate at adjacent time points, construct a joint transfer matrix and interpolate between adjacent time points to obtain gene expression data over continuous time. S2. Based on continuous-time gene expression data, an ordinary differential equation model incorporating transcription factor expression levels is used to output global transcription factor-gene regulatory weights as the first prior. S3. Based on single-cell chromatin openness data, construct a multi-omics association matrix and generate a transcription factor-gene regulation prior matrix. Then, fuse the prior matrix with the first prior to obtain the enhanced prior. S4. Determine the graph topology of the graph neural network based on the enhanced prior, and output the transcription factor-gene regulatory weight matrix and cell phenotype prediction value at each intermediate time point, with transcription factors and genes as nodes.

[0006] Preferably, S1 further includes: The cost matrix is ​​defined based on the gene expression distance in the scRNA-seq data of adjacent time points. The optimal transfer problem is solved to obtain the optimal transfer matrix. The optimal transfer matrix is ​​then fused with the local transfer matrix derived from the RNA rate to obtain the joint transfer matrix. Based on the joint transfer matrix, the cell states at multiple intermediate time points are interpolated between adjacent time points to obtain continuous-time gene expression, splicing, and unspliced ​​data.

[0007] Preferably, S2 further includes: Based on the cell2fate ordinary differential equation model framework, the modular transcription rate formula is modified by introducing the contribution of transcription factor expression level to the modular transcription rate. According to the module activation time and module time weight, the total regulatory weight is output as the global transcription factor-gene regulatory weight, i.e., the first prior.

[0008] Preferably, S3 further includes: Open chromatin peaks were identified using scATAC-seq data, and the openness score of each peak was calculated. A time-specific peak-gene association matrix was constructed, and a transcription factor-peak binding matrix was generated by scanning transcription factor binding motifs in the peak sequences. Based on the transcription factor-peak binding matrix and peak openness scores, a transcription factor-gene regulation prior matrix was generated. Then, the transcription factor-gene regulation prior matrix was fused with the first prior using a Hadamard product to obtain the enhanced prior.

[0009] Preferably, S4 further includes: Using transcription factors and genes as nodes, node features are integrated with information related to gene expression levels, transcription factor expression levels, and chromatin accessibility. The initial graph topology is determined by thresholding based on enhancement priors. The GraphSAGE layer is used to aggregate neighbor information, and the time-resolved transcription factor-gene regulatory weight matrix is ​​output through edge MLP and sigmoid function.

[0010] Preferably, the graph neural network in S4 is a GraphSAGE+MLP graph neural network, and its loss function simultaneously constrains: the transcription rate predicted by the model is consistent with the RNA rate in S1, phenotypic prediction error, L1 sparse regularization, and the model weights approach the enhanced prior.

[0011] Preferably, in S1, the single-cell transcriptome data includes expression levels, splicing levels, and unspliced ​​levels; and the chromatin accessibility data is not interpolated and is aligned to the most recent time point.

[0012] Preferably, in step S3, when there is no scATAC-seq data, the transcription factor-gene regulation prior matrix is ​​an identity matrix, and the enhanced prior is equivalent to the first prior.

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention enables precise reconstruction of the continuous evolution of gene regulatory networks from discrete single-cell multi-omics data, effectively avoiding directional biases at key differentiation nodes. By fusing transcriptomic and chromatin accessibility information, the reliability of regulatory priors is significantly improved, allowing for accurate differentiation of complex regulatory scenarios. The end-to-end graph neural network architecture mutually reinforces regulatory inference and phenotypic prediction, and error signals can be backpropagated to optimize regulatory relationships, thereby obtaining a high-resolution, high-precision dynamic regulatory weight matrix. This method remains robust even with sparse data, adaptively balancing priors and exploration at key bifurcation states, significantly reducing false negative and false positive biases, and providing reliable computational support for analyzing dynamic processes such as development and disease. Attached Figure Description

[0014] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 This is a schematic diagram of the model structure according to an embodiment of the present invention; Figure 2 This is a visualization of UMAP in the cell type-specific regulatory network analysis of mouse erythroid cells in an embodiment of the present invention. Figure 3 This is a diagram of the Sp1 regulator network in the cell type-specific regulatory network analysis of mouse erythroid cells in an embodiment of the present invention. Figure 4 This is a schematic diagram of the top 20 differentially expressed genes in five cell types in the cell type-specific regulatory network analysis of mouse erythroid cells in an embodiment of the present invention. Figure 5 This is a distribution diagram of the expression of six genes regulated by Sp1 in the cell type-specific regulatory network analysis of mouse erythroid cells in an embodiment of the present invention. Figure 6 This is a heatmap showing the average activity of eight gene expression modules in five cell types during signal flow analysis of mouse erythroid cells in an embodiment of the present invention. Figure 7 This is a heatmap showing the average activity of eight GEMs during developmental stages in the signal flow analysis of mouse erythroid cells according to an embodiment of the present invention. Figure 8 This is a volcano plot showing the differential expression of target genes E7.0 vs E8.5 in mouse erythroid cells during signal flow analysis in an embodiment of the present invention. Figure 9 This is a signal flow network diagram at time points E7.0, E8.0, and E8.5 in the signal flow analysis of mouse erythroid cells according to an embodiment of the present invention. Figure 10 This is a Venn diagram of the signal flow regulation edges at E7.0 and E8.5 in the signal flow analysis of mouse erythroid cells in an embodiment of the present invention. Figure 11 Venn diagrams of the target gene set in the relationship between GEM and target genes at E7.0 and E8.5 in the signal flow analysis of mouse erythroid cells in this embodiment of the invention. Figure 12 This is a time-series dynamic graph showing the number of three types of signal flow edges during the developmental stage in the signal flow analysis of mouse erythroid cells according to an embodiment of the present invention; Figure 13This is a UMAP dimensionality reduction visualization of the MouseErythroid development dataset in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention (colored according to the Louvain algorithm clustering results). Figure 14 This is a bubble diagram showing the degree centrality of the first 20 TFs in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention. Figure 15 This is a degree centrality ranking diagram of the Top 30 TFs in the Blood progenitors 1 and Erythroid 3 cell types in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention. Figure 16 The diagram shows the cell identity displacement vector field and randomized control vector field after Tead1 gene knockout simulation in the cell perturbation analysis of mouse erythroid cells in this embodiment of the invention. Figure 17 This is a spatial distribution diagram on UMAP of the inner product score of the KO displacement vector and the natural differentiation direction vector of each cell in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention. Figure 18 This is a spatial distribution of cell pseudotime and natural differentiation trajectory vector field on UMAP in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention; Figure 19 This is a high-resolution visualization on UMAP of the inner product of the KO displacement vector and the natural differentiation direction in the cell perturbation analysis of mouse erythroid cells in an embodiment of the present invention. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0017] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0018] Example 1 This embodiment provides a dynamic gene regulatory network inference method based on optimal transport and multi-omics prior fusion, the steps of which include: S1. Based on single-cell transcriptome data and RNA rate at adjacent time points, construct a joint transfer matrix and interpolate between adjacent time points to obtain gene expression data over continuous time.

[0019] First, OT transfer matrices are constructed from scRNA-seq data at adjacent time points, including expression levels, splicing levels, and unspliced ​​levels. These are then combined with local transfer matrices derived from RNA rates to generate joint transfer matrices. This allows for the interpolation of cell states at multiple intermediate time points between adjacent time points, resulting in continuous-time gene expression, splicing, and unspliced ​​data. ATAC data are not interpolated and are aligned according to the most recent time point.

[0020] Specifically, take two adjacent time points. and The scRNA-seq data includes expression data. X splicing data S Uncut data U ,in n for The number of cells at any given time. m for The number of cells at any given time. g Number of genes, and on X , S , U Z-score standardization preprocessing was performed separately, followed by OT computation. Every cell at any moment The correspondence between each cell at each time step. This embodiment first defines the cost matrix. , where each element Ci,j express cells i and cells The gene expression distance between them, and the cost matrix is ​​shown in Equation (1): (1) (2) Next, the optimal transmission problem is solved, that is, the transmission matrix is ​​obtained through formula (2). Then, scVelo was used to calculate the RNA rate and construct a local transition matrix. By combining the two, we obtain the joint transition matrix: (3) in, Balance OT and Normalization make sure .

[0021] (4) (5) Finally, this embodiment considers adjacent time points. t and Interpolate the cell states at multiple intermediate time points. However, due to a cell i It may evolve into multiple cells. j (Probability distribution), take the expected expression vector of cell i at time s as shown in formula (5), and finally obtain gene expression, splicing and unspliced ​​data for continuous time; ATAC data are aligned according to the most recent time point and no interpolation is performed.

[0022] S2. Based on continuous-time gene expression data, an ordinary differential equation model incorporating transcription factor expression levels is used to output global transcription factor-gene regulatory weights as the first prior.

[0023] Based on the cell2fate-based ODE model framework, the contribution of TF expression level to modular regulatory rate is introduced, and the global TF-gene regulatory weights are output as priors for subsequent graph neural networks.

[0024] Specifically, based on unspecified u splicing s mRNA expression levels were determined by fitting rate parameters using an ODE model for each cell. c and genes g This embodiment has the following ODE model: (6) (7) (8) in, The rate at which module m reaches the target transcription rate in state i. Representative module m In state i Lower gene g The target transcription rate, For module m In state i Convergence rate under (ON / OFF) conditions.

[0025] Then introduce To modify the modular transcription rate formula (8) in the known model, we obtain formula (9).

[0026] (9) in, For module m Transcription factors f On genes g The regulatory weight, The expression level of transcription factor f, then according to: (10) (11) (12) Formula (10) represents the module activation time, formula (11) represents the module time weight, and formula (12) represents the overall regulatory weight. The final output is the RNA rate, along with the TF-gene weights as a global prior. .

[0027] S3. Based on single-cell chromatin accessibility data, construct a multi-omics association matrix and generate a transcription factor-gene regulation prior matrix. Then, fuse the prior matrix with the first prior to obtain the enhanced prior.

[0028] Time-specific peak-gene association matrices and TF-peak binding matrices were constructed using scATAC-seq data. Combined with peak openness scores, a TF-gene regulation prior matrix was generated. The enhanced prior is then fused with the Hadamard product of the output prior of S2 to obtain the enhanced prior. .

[0029] Specifically, scATAC-seq data was used to identify open chromatin peaks, annotate the genomic regions of the peaks, calculate the openness score for each peak, and associate the peaks with genes. Then, a time-specific peak-gene association matrix was constructed. Scan the TF-binding motifs in the Peak sequence to generate the TF-peak binding matrix. Combined with peak openness score Generate the ATAC prior matrix: (13) Will The enhanced prior is obtained by fusing the prior weights of S2 with the Hadamard product: (14) When there is no scATAC-seq data .

[0030] S4. Determine the graph topology of the graph neural network based on the enhanced prior, and output the transcription factor-gene regulatory weight matrix and cell phenotype prediction value at each intermediate time point, with transcription factors and genes as nodes.

[0031] A GraphSAGE+MLP graph neural network is constructed using TFs and genes as nodes. Node features are fused from gene expression, TF expression, ATAC-related peak openness, and binding information. The graph topology is determined by a threshold based on fusion priors. The loss function simultaneously constrains: ( a The model predicted a transcription rate consistent with the RNA rate in S1; bPhenotypic prediction error; c L1 sparse regularization; d The model weights approximate the fusion priors. The final output is the TF-gene regulatory weight matrix and cell phenotype prediction values ​​for each intermediate time point.

[0032] Specifically, using transcription factors and genes as nodes, node features are fused to determine expression levels and ATAC-related information; based on fusion priors... The initial graph topology is determined by a threshold; neighbor information is aggregated using a GraphSAGE layer; and time-resolved TF-genes are used to regulate weights through edge MLP and sigmoid function outputs. Further computation of cell-level embeddings and prediction of phenotypes. The model's loss function simultaneously constrains transcription rate consistency, phenotypic prediction error, L1 sparse regularization, and prior constraints. The loss function is as follows: (15) The first term is the mean squared error (MSE), which ensures the transcription rate predicted by the model. RNA rate with S1 Consistency; if the second term is a categorical phenotype, such as cell type or functional state, then cross-entropy loss is used; if it is a continuous phenotype, such as the degree of differentiation of continuous values ​​between 0 and 1, then mean squared error is used; the third term is L1 regularization, which promotes the sparsity of GRN; the fourth term is prior constraint, which ensures that the model weights are close to the prior weights. It is the regularization coefficient, which is adjusted through cross-validation.

[0033] The implementation of this embodiment is based on, as follows Figure 1 The network model shown includes: an optimal transport OT interpolation module, an RNA rate-driven TF-gene regulation prior module, a multi-omics prior fusion module, and a time-resolved GRN inference and phenotypic prediction module. These modules work together to achieve the steps described in the above embodiments.

[0034] The following will describe in detail how the present invention solves the technical problems in practical work, using mouse erythroid cells as an example.

[0035] (1) Analysis of cell type-specific regulatory network of mouse erythroid cells Traditional transcriptional regulatory network analysis, often based on batch cell sequencing data, struggles to elucidate the unique regulatory logic of different cell types within tissues, frequently obscuring cell type-specific transcription factor activity and their target gene modules. This embodiment introduces a cell type-specific regulatory network analysis strategy. First, it clearly defines the target cell population, then integrates transcription factor binding motifs and co-expression networks to construct a highly reliable cell type-specific regulatory subset. By comparing the regulatory networks of two cell types, this embodiment aims to screen for core transcription factors that determine cell identity and functional specialization.

[0036] Based on the regulatory edges obtained from the model in this embodiment, true regulatory edges with scores greater than 0.01 were selected. Erythroid2 had the most regulatory edges, with 9706. This embodiment considers that the higher the degree of a transcription factor, the more important it is. Sp1 regulates 523 target genes and can be considered a core transcription factor. Figure 2 Visualization of five cell type populations based on UMAP, including blood progenitors (Blood progenitors 1 / 2) and erythrocyte differentiation populations (Erythroid 1 / 2 / 3). Figure 3 This represents the regulatory subnetwork of Sp1. Figure 4 This represents the top 20 differentially expressed genes for these five cell types. This example further examined genes regulated by Sp1 and found that six genes, including Hba-x, Mif, and Hba-a2, are also marker genes in Erythroid3 cells. This example then examined the expression of these six genes, from... Figure 5 They were observed to be highly expressed in Erythroid2 cells and lowly expressed in other cell types. Hba-x, Mif, and Hba-a2 play key regulatory roles in the Erythroid2 gene regulatory network.

[0037] (2) Signal flow analysis of mouse erythroid cells Intercellular communication drives directional information flow triggered by gene expression modules within cells, which in turn triggers the outflow of other signals. Gene expression modules (GEMs), driven by the regulatory relationship between transcription factors and target genes, can also form directional information flow within the cell. This embodiment utilizes NMF dimensionality reduction to obtain gene modules (GEMs), and then, based on OLS regression combined with biological prior constraints, infers the causal directed network of "transcription factor → gene module → target gene," thereby revealing the regulatory information flow pathways during development.

[0038] Based on the control edges obtained from the model in this embodiment, real control edges with scores greater than 0.01 are selected for signal flow analysis. Figure 6 It can be seen that, on average, different cell types activate different GEM combinations, revealing that cell identity is determined by the synergistic action of specific modules. Figure 7 On average, during the developmental stage, the early stage (E7.0-E7.5) is dominated by a specific GEM, the late stage (E8.0-E8.5) sees the emergence of new module activation patterns, and the middle stage is a transitional state in which module activity changes. Figure 8 Volcano plots revealed differential expression of target genes at two developmental poles, identifying representative differentially expressed genes Pa2g4, Clint1, Helz, and Atg16l1. These differentially expressed genes reflect the molecular programming transition from early to late erythrocyte development. Pa2g4 was significantly upregulated (p=9.5e-14, log2FC≈1.2); Clint1 was significantly upregulated (p=0.0028, log2FC≈0.96); Helz was marginally significant (p=0.062, log2FC≈2.4); and Atg16l1 showed a downregulation trend (p=0.095, log2FC≈-0.93). Figure 9 The temporal dynamics of the signal flow network are shown. The scale of the signal flow expands significantly with the developmental process. The network size is smallest in E7.0 (early stage), with only 8 inflow edges, 624 outflow edges, and 10 GEM-GEM interactions. The regulatory program is relatively simple. The number of cells is the largest in E8.0 (mid-stage peak) (2,811), with 333 inflow edges and 1,633 outflow edges. The complexity of the regulatory network reaches its peak. E8.5 (late stage) has the largest network size, with 367 inflow edges, 1,931 outflow edges, and 24 GEM-GEM interactions. The coordination between modules is the most complex. Figure 10 Venn diagram analysis revealed a conserved signal flow edge between E7.0 and E8.5, as well as a large number of period-specific regulatory events. The conserved edge may represent the core essential program of erythrocyte development, while the specific edge reflects the regulatory needs of developmental stages. Figure 11 The Venn analysis of GEM target genes also showed that some target genes were regulated in both phases, with a new set of target genes appearing in the later stage, suggesting a functional shift. Figure 12 The time-series dynamics of the number of three types of signal flow edges during the development period.

[0039] (3) Cell perturbation analysis of mouse erythroid cells To predict the function of key transcription factors (TFs) during development, this embodiment simulates the insilicoTead1 gene knockdown (KD) perturbation based on a constructed temporal GRN. In the GRN, Tead1 regulates 174 target genes. This embodiment simulates the KD effect by assuming a 50% reduction in target gene expression and projects the expression changes into a low-dimensional UMAP space to calculate the displacement vector for each cell. Cell population annotation and TF centrality analysis are also performed.

[0040] Figure 13The Louvain clustering annotation results for 9,815 cells clearly distinguish the continuous differentiation trajectories of hematopoietic stem cell progenitors (Bloodprogenitors1 / 2) and erythrocyte lineages (Erythroid1 / 2 / 3). Figure 14 The bubble diagram shows the weighted out-degree and degree-centrality distribution of TFs, with TFs such as Trp53, Sp1, and Ets1 located in regions of high regulatory capacity, suggesting their pivotal role in the erythrocyte development network. Figure 15 Further comparison of the ranking of the top 30 core TFs in hematopoietic stem cell progenitors (Bloodprogenitors1) and terminal erythroid cells (Erythroid3) revealed that Trp53 ranked first in both populations, while Nfe2l2, Stat3, and others were significantly enriched in the erythroid stage, reflecting the lineage-specific dynamic changes in TF activity. Figure 16 The cell displacement vector field after Tead1KD (left) and the random control (right) are shown. Compared with random perturbation, Tead1KD induced significant directional displacement in the erythrocyte region, with a mean displacement amplitude of 0.0267, indicating that Tead1 loss does indeed alter the transcriptomic state of cells. To determine whether Tead1KD promotes or inhibits erythrocyte differentiation, this embodiment calculated the normalized inner product of the KO displacement vector and the natural differentiation gradient. Figure 17 The results showed that the inner volume value of most cells was close to 0, but local negative values ​​(blue) appeared in the erythrocyte region (Erythroid2 / 3), suggesting that Tead1KD may delay or arrest terminal differentiation of erythrocytes. Figure 18 The pseudo-temporal distribution and natural differentiation direction vector field were further demonstrated, verifying the differentiation axis from progenitor cells to terminal erythrocytes. Figure 19 Mapping the inner product to the perturbation vector field reveals that in late-differentiation cell populations, the KO direction shows an inverse trend compared to the natural differentiation direction, further supporting the hypothesis that Tead1 is a positive regulator of erythrocyte differentiation. In summary, insilico perturbation analysis based on time-series GRNs demonstrates that Tead1 plays a crucial positive role in terminal erythrocyte differentiation by regulating 174 target genes; its deletion can lead to cells deviating from the normal differentiation trajectory, suggesting that Tead1 is a key regulator in erythrocyte development.

[0041] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion, characterized in that, Includes the following steps: S1. Based on single-cell transcriptome data and RNA rate at adjacent time points, construct a joint transfer matrix and interpolate between adjacent time points to obtain gene expression data over continuous time. S2. Based on continuous-time gene expression data, an ordinary differential equation model incorporating transcription factor expression levels is used to output global transcription factor-gene regulatory weights as the first prior. S3. Based on single-cell chromatin openness data, construct a multi-omics association matrix and generate a transcription factor-gene regulation prior matrix. Then, fuse the prior matrix with the first prior to obtain the enhanced prior. S4. Determine the graph topology of the graph neural network based on the enhanced prior, and output the transcription factor-gene regulatory weight matrix and cell phenotype prediction value at each intermediate time point, with transcription factors and genes as nodes.

2. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion as described in claim 1, characterized in that, S1 further includes: The cost matrix is ​​defined based on the gene expression distance in the scRNA-seq data of adjacent time points. The optimal transfer problem is solved to obtain the optimal transfer matrix. The optimal transfer matrix is ​​then fused with the local transfer matrix derived from the RNA rate to obtain the joint transfer matrix. Based on the joint transfer matrix, the cell states at multiple intermediate time points are interpolated between adjacent time points to obtain continuous-time gene expression, splicing, and unspliced ​​data.

3. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion according to claim 1, characterized in that, S2 further includes: Based on the cell2fate ordinary differential equation model framework, the modular transcription rate formula is modified by introducing the contribution of transcription factor expression level to the modular transcription rate. According to the module activation time and module time weight, the total regulatory weight is output as the global transcription factor-gene regulatory weight, i.e., the first prior.

4. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion according to claim 1, characterized in that, S3 further includes: Open chromatin peaks were identified using scATAC-seq data, and the openness score of each peak was calculated. A time-specific peak-gene association matrix was constructed, and a transcription factor-peak binding matrix was generated by scanning transcription factor binding motifs in the peak sequences. Based on the transcription factor-peak binding matrix and peak openness scores, a transcription factor-gene regulation prior matrix was generated. Then, the transcription factor-gene regulation prior matrix was fused with the first prior using a Hadamard product to obtain the enhanced prior.

5. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion according to claim 1, characterized in that, S4 further includes: Using transcription factors and genes as nodes, node features are integrated with information related to gene expression levels, transcription factor expression levels, and chromatin accessibility. The initial graph topology is determined by thresholding based on enhancement priors. The GraphSAGE layer is used to aggregate neighbor information, and the time-resolved transcription factor-gene regulatory weight matrix is ​​output through edge MLP and sigmoid function.

6. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion according to claim 1, characterized in that, The graph neural network in S4 is a GraphSAGE+MLP graph neural network, and its loss function simultaneously constrains: the transcription rate predicted by the model is consistent with the RNA rate in S1, phenotypic prediction error, L1 sparse regularization, and the model weights approach the enhanced prior.

7. The method for inferring dynamic gene regulatory networks based on optimal transport and multi-omics prior fusion according to claim 1, characterized in that, In S1, the single-cell transcriptome data includes expression levels, splicing levels, and unspliced ​​levels; and the chromatin accessibility data is not interpolated and is aligned to the most recent time point.

8. The method according to claim 1, characterized in that, In step S3, when there is no scATAC-seq data, the transcription factor-gene regulation prior matrix is ​​an identity matrix, and the enhanced prior is equivalent to the first prior.

Images