A gene regulation inference method guided by topological data analysis for gene network embedding

Abstract

This invention discloses a gene regulation inference method guided by topological data analysis and gene network embedding. It combines TDA and GNN to enhance the inference capability of gene regulation networks. By capturing the topological structure of the gene regulation network graph through TDA features, the model's ability to model gene expression is enhanced. The TDA features and GAT embedding representations are effectively integrated through gating fusion. This fusion mechanism enables the model to adaptively adjust node embeddings based on global topological characteristics, which not only improves the accuracy of gene interaction representation but may also enhance the accuracy of regulatory relationship prediction. The traditional GAT architecture is extended through a four-layer graph attention mechanism. Each layer uses residual connections to alleviate the gradient vanishing problem and improve training stability. In addition, independent multilayer perceptron branches are designed for transcription factors and target gene embeddings. This deep architecture can achieve more expressive feature transformations and capture subtle patterns in gene regulation networks.

Description

Technical Field

[0001] This invention relates to the technical field of bioinformatics, and in particular to a method for inferring gene regulation by guiding gene network embedding through topological data analysis. Background Technology

[0002] Gene regulatory networks, or GRNs, characterize the complex regulatory relationships between transcription factors and their target genes. The advent of single-cell RNA sequencing technology has revolutionized the study of cellular heterogeneity, making it possible to perform high-resolution gene expression analysis at the single-cell level. Unlike batch RNA sequencing, which averages population signals and thus masks cell-specific variations, scRNA-seq data can capture unprecedented details of gene interactions, providing a unique opportunity to infer cell type-specific GRNs. However, this fine-grained data also brings significant challenges, including high sparsity due to loss events, measurement noise, cellular heterogeneity, and the inherent nonlinearity of regulatory interactions. These factors make accurate GRN inference more complex.

[0003] Chinese patent application CN202211571759.6 discloses a method for inferring gene regulatory networks based on a linear mixture model. This method takes single-cell gene expression data as input and represents the data using mean, random effects, and noise. The random effects and noise are both matrix random variables that follow a normal matrix distribution. Then, the row covariance matrix is ​​obtained through parameter estimation to represent the correlation between genes. However, this method still has shortcomings. It only considers steady-state data and does not simultaneously consider the temporal information contained in the expression data. This results in incomplete single-cell gene expression data used before inference, failing to provide dynamic regulatory information about genes and affecting the accuracy of the inferred gene regulatory network. Furthermore, by focusing only on the data, it to some extent ignores the inherent topological characteristics of the gene regulatory network, leading to a lack of clear biological explanation for its structure.

[0004] Chinese patent application CN201911012024.8 discloses a method for constructing a gene regulatory network based on scRNA-seq and dynamic time warping, including the following steps: Step 1, single-cell pseudo-time analysis based on single-cell data; Step 2, dual DTW distance analysis based on positive and negative regulatory modes; Step 3, cluster analysis and GRN module screening based on DTW distance; Step 4, construction of a gene similarity matrix based on DTW distance; Step 5, calculation of the gene-gene network connectivity threshold C through random permutation test; Step 6, gene functional analysis of GRN modules; Step 7, visualization of the pGRN network to obtain the GRN network diagram, which can detect time-regulated... This method addresses both time-lag and time-scaling regulatory relationships. By applying the DTW algorithm, it can handle similar expression map searches in time-scaling scenarios, expanding the search scope for regulatory relationships over time. However, this method still has shortcomings. It only considers the information contained in time-series data, ignoring the gene regulatory relationships in steady-state data. The single-cell gene expression data used before inference is not complete enough and is easily affected by noise, making it impossible to guarantee the stability of the network. Although clustering is performed to obtain corresponding modules before inferring the gene regulatory network, it only uses simple similarity as the basis for whether there are relationships between genes within the modules, resulting in an inaccurate and unreliable gene regulatory network.

[0005] Most existing gene regulatory network inference methods rely on various machine learning and deep learning-based approaches to infer single-cell gene regulatory networks. While these methods have demonstrated their respective advantages, they often focus on information aggregation strategies for nodes and edges within graph data, neglecting the more complex topological features in high-dimensional graph data and failing to extract relevant information more comprehensively from the graph structure. Summary of the Invention

[0006] The purpose of this invention is to solve the problems in the prior art and propose a gene regulation inference method guided by topological data analysis and gene network embedding. This method can combine a global topological perspective and graph local representation, capture the topological structure of the gene regulation network graph through TDA features, thereby enhancing the model's ability to model gene expression. By using the enhanced gene expression embedding, the identification of interactions between genes can be improved.

[0007] To achieve the above objectives, this invention proposes a method for inferring gene regulation through gene network embedding guided by topological data analysis, comprising the following steps:

[0008] S1. Topological Data Analysis Feature Calculation: Calculate TDA features from gene expression data and prior knowledge to capture the global topological characteristics of the gene regulatory graph. The persistence value of the Betti number reflects the change of topological features during the graph filtering process. The persistence value is obtained by using a filtering function. The persistence process of the graph is constructed using persistent graphs and persistent barcodes. The TDA feature embedding is obtained using persistent barcodes.

[0009] S2. Feature encoder based on graph attention network: A multi-layer GAT layer is used as the feature encoder of TDAGENE. Each GAT layer introduces a self-attention mechanism in the information propagation step, and assigns importance weights to the neighbors in the graph through the nodes.

[0010] S3, Fusion Layer: A gated fusion mechanism is used to integrate GAT embeddings with TDA features to generate enhanced representations;

[0011] S4, Branch-specific multilayer perceptron: The fused embedding is processed through independent transcription factor branches and target gene branches;

[0012] S5, Decoder: The final regulatory value is calculated using an MLP decoder. Based on the transcription factor embedding and the target gene embedding, the embeddings of the two branches are spliced ​​together and processed through a three-layer neural network. Each layer contains linear transformation, layer normalization and LeakyReLU to decode and predict the regulatory links between gene pairs.

[0013] S6. Training Optimization: Use binary cross-entropy as the loss function to measure the error between the predicted result and the true label, and use its gradient to update the network weights to drive continuous optimization of model performance.

[0014] Preferably, in step S1, a graph structure is constructed using prior knowledge, and a filtering function is used to identify global topological features of the graph from the graph structure, including connected components and one-dimensional holes. Persistent values ​​reflect changes in global topological features, and persistent value graphs and persistent value barcodes are used to obtain the embedding of global topological features.

[0015] Preferably, in step S1, the topological features of the graph are extracted and represented through the graph filtering process, modeled using a neural network, and the graph learning task is trained using a filtering function to optimize the filtering process. The filtering function is based on nodes and edges, and TDA feature embeddings are obtained by calculating persistent homology for a given graph, thereby extracting global multi-scale information of the graph dataset.

[0016] Preferably, in step S2, the feature encoder takes the sparse adjacency matrix and gene expression matrix representing prior regulatory relationships as input, splices the normalized node degree as an additional feature, and splices the input dimension with the additional feature dimension as the effective input dimension.

[0017] Preferably, there are four GAT layers in step S2. The output of each multi-head attention module is aggregated based on the mean or splicing. Each layer introduces residual connections with linear projection to generate GAT embeddings that capture local and contextual gene interactions.

[0018] Preferably, in step S3, the TDA vector of each node is repeated and adjusted to match the dimension of the GAT embedding, forming an expanded topological feature matrix. Then, the GAT embedding is concatenated with the topological features, and a weight gate vector is generated through a Sigmoid gate unit. The embedding vector is adjusted through this gate vector, and the fused output is further refined by the projection MLP so that the embedding contains topological information.

[0019] Preferably, each branch in step S4 contains three linear layers, and each module sequentially includes linear transformation, layer normalization, LeakyReLU activation, and Dropout regularization. The final embedding vector of the transcription factor branch is output through a projection layer, and the target gene branch is processed using a symmetrical structure. Specific feature learning is performed for different regulatory roles of transcription factors and target genes, and LeakyReLU is used as the activation function of the output layer.

[0020] Preferably, in step S6, a full-batch gradient descent strategy is adopted, a reparameterization technique is introduced for gradient estimation, the Adam optimizer is used, the learning rate is set to 0.005, the weight decay coefficient is 0.0005, and all key hyperparameters are optimized by grid search.

[0021] The beneficial effects of this invention are as follows: This invention combines a global topological perspective with local graph representation, integrating topological data analysis into a graph-based neural network framework to capture local and global structural information in gene regulatory networks and predict regulatory links. Although graph topological features cannot be directly modeled, they can be extracted using the Union-Find algorithm through persistent cohomology computation based on cosine similarity filtering values ​​and fused with graph embeddings to improve GRN inference. This can be described by incorporating TDA features into a multi-layer graph attention network with residual connections, constructing a topological mechanism framework, TDAGENE, which provides enhanced accuracy GRN inference through topological insights. TDAGENE takes a sparse adjacency matrix A representing prior regulatory relationships and a gene expression matrix X as inputs. It utilizes a multi-head GAT layer as an encoder to generate node embeddings from gene features. These embeddings are then fused with global TDA features using a gated fusion mechanism. The fused representation is passed through a two-branch MLP to produce specialized embeddings. These learned embeddings are then input to a decoder layer, which uses an MLP to decode and predict regulatory links between gene pairs. In this study, TDAGENE employs residual projection and layer normalization to stabilize deep GAT propagation, ensuring robust feature learning on sparse biological graphs. It captures the topological structure of gene regulatory network graphs through TDA features, thereby enhancing the model's ability to model gene expression. Enhanced gene expression embeddings improve the identification of interactions between genes. Combining TDA and GNNs enhances the inference ability of gene regulatory networks. TDA features capture the topological structure of gene regulatory network graphs, further enhancing the model's ability to model gene expression. Gated fusion effectively integrates TDA features and GAT embeddings, enabling the model to adaptively adjust node embeddings based on global topological characteristics. This not only improves the accuracy of gene interaction representation but may also enhance the accuracy of regulatory relationship prediction. A four-layer graph attention mechanism extends the traditional GAT architecture, with residual connections in each layer to alleviate the gradient vanishing problem and improve training stability. Furthermore, independent multilayer perceptron branches are designed for transcription factors and target gene embeddings. This deep architecture enables more expressive feature transformations and can capture subtle patterns in gene regulatory networks.

[0022] The features and advantages of the present invention will be described in detail through embodiments and in conjunction with the accompanying drawings. Attached Figure Description

[0023] Figure 1 This is a TDAGENE framework diagram of a gene regulation inference method guided by topological data analysis and gene network embedding according to the present invention;

[0024] Figure 2This is a schematic diagram illustrating the use of persistent graphs to retain the coherence information of a given graph through a graph filtering process in the TDA feature calculation framework of a gene regulation inference method guided by topological data analysis and gene network embedding according to the present invention.

[0025] Figure 3 This is a schematic diagram of the TDAGENE framework filtering learning process of a gene regulation inference method guided by topological data analysis and gene network embedding according to the present invention.

[0026] Figure 4 This is an AUROC score graph of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in this invention on a non-specific ChIP-seq dataset.

[0027] Figure 5 This is an AUROC score map of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in this invention, on a cell type-specific ChIP-seq dataset.

[0028] Figure 6 This is the AUROC score map of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in the present invention on the STRING protein interaction dataset.

[0029] Figure 7 This is an AUPRC score map of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in this invention, on a non-specific ChIP-seq dataset.

[0030] Figure 8 This is an AUPRC score map of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in this invention, on a cell type-specific ChIP-seq dataset.

[0031] Figure 9 This is the AUPRC score map of 10 models of a gene regulation inference method guided by topological data analysis and gene network embedding in the present invention on the STRING protein interaction dataset.

[0032] Figure 10 This is a UMAP dimensionality reduction visualization result of the gene embedding vector output by TDAGENE on the human embryonic stem cell hESC dataset, which is a gene regulation inference method guided by topological data analysis for gene network embedding according to the present invention.

[0033] Figure 11 This is a dynamic change diagram of gene regulatory loops in a time series identified by the TDAGENE method of the present invention, which is a gene regulation inference method for gene network embedding guided by topological data analysis.

[0034] Figure 12 This is a graph showing the changes in gene expression levels of three core transcription factors, NANOG, SOX2, and POU5F1, in human embryonic stem cells, based on a gene regulation inference method guided by topological data analysis and gene network embedding according to the present invention.

[0035] Figure 13 This is a visualization of a gene regulation inference method based on topological data analysis-guided gene network embedding, based on the present invention, using TDAGENE to predict real cell type networks.

[0036] Figure 14 This is a graph showing the degree distribution and main topological properties of a gene regulation network constructed on the hESC dataset by TDAGENE, a gene regulation inference method for gene network embedding guided by topological data analysis according to the present invention.

[0037] Figure 15 This is a graph showing the number of times transcription factor TF is regulated in the real prior regulatory network of a gene regulation inference method that guides gene network embedding by topological data analysis according to the present invention.

[0038] Figure 16 This is a graph showing the impact of four hyperparameter values ​​on model performance of a gene regulation inference method based on topological data analysis-guided gene network embedding in this invention.

[0039] Figure 17 This invention presents four variant models of a gene regulation inference method guided by topological data analysis and gene network embedding, and AUROC value plots of TDAGENE on datasets of seven different cell types.

[0040] Figure 18 This invention presents four variant models of a gene regulation inference method guided by topological data analysis and gene network embedding, and AUPRC value plots of TDAGENE on datasets of seven different cell types. Detailed Implementation

[0041] 1.1 Problem Definition

[0042] Given a gene expression profile of an scRNA-seq dataset with N cells and M genes. The primary goal is to infer the potential regulatory relationships between transcription factors and target genes. By utilizing this prior knowledge, the task of inferring gene regulatory networks is transformed into a network connectivity prediction problem. This prior knowledge is represented by G(V, E) in the gene regulation graph, where V represents the set of M gene nodes, E represents the set of observed edge regulatory relationships, and G is related to the adjacency matrix. Associated, where only if (vi,vj)∈E ,otherwise In the gene expression spectrum X, each column represents the node features of the corresponding gene. The goal of this invention is to predict the missing edges in the gene regulatory network graph G.

[0043] 1.2 Model Framework

[0044] This invention aims to integrate topological data analysis into a graph-based neural network framework to capture local and global structural information in gene regulatory networks and predict regulatory links. While graph topological features cannot be directly modeled, they can be extracted using the Union-Find algorithm through persistent cohomology computation based on cosine similarity filtering values ​​and fused with graph embeddings to improve GRN inference. This can be described by incorporating TDA features into a multi-layer graph attention network with residual connections, such as... Figure 1 As shown, to this end, the present invention constructs a topological mechanism framework, namely TDAGENE, which provides GRN inference with enhanced accuracy through topological insights. TDAGENE takes a sparse adjacency matrix A representing prior regulatory relationships and a gene expression matrix X as inputs. It uses a multi-head GAT layer as an encoder to generate node embeddings from gene features. These embeddings are then fused with global TDA features computed from the graph using a gated fusion mechanism. The fusion representation is passed through a two-branch MLP, namely a transcription factor branch and a target gene branch, to generate specialized embeddings. These learned embeddings are then input to a decoder layer, which uses an MLP to decode and predict the regulatory links between gene pairs. Between multiple GAT layers, TDAGENE uses residual projection and layer normalization to stabilize deep GAT propagation, ensuring robust feature learning on sparse biological graphs.

[0045] Figure 1 The TDAGENE framework is presented. (a) TDAGENE extracts GAT and TDA features from gene expression data and prior knowledge. By fusing the two features, more comprehensive multi-scale information in the graph is obtained. (b) Graph attention with residual connections captures the interaction between transcription factors and genes using the attention mechanism. Combined with residual connections, the stability of training is ensured. (c) GAT embeddings and TDA features are fused through gating fusion. The fused features are then processed by a branch-specific multilayer perceptron for specialized feature learning. Finally, the obtained features are concatenated and processed by an MLP decoder to obtain the regulatory relationship score between genes.

[0046] 1.3 TDA Feature Calculation

[0047] This invention captures the global topological properties of gene regulatory graphs by calculating TDA features. From a topological perspective, graphs contain topological features of 0-dimensional connected components and 1-dimensional cycles. This invention can effectively extract and represent the topological features of graphs by using neural networks for modeling through a graph filtering process, such as... Figure 2 As shown, this invention utilizes a filtering function to train a graph learning task to optimize the filtering process. The filtering function is based on nodes and edges, and the computation process of this invention begins with node embedding. and Cosine similarity is used to construct graph edge filtering, which is defined as:

[0048] (1);

[0049] in As an edge The filter value, considering only the upper triangular edges in the sparse adjacency matrix index, is used to construct the persistent homology filtering process of the graph, representing the "distance" between nodes i and j based on cosine similarity.

[0050] This invention arranges the edges in ascending order of their filter values: ,satisfy Subsequently, continuous homology is applied through a union-find structure to extract persistent pairs of dimensions 0 and 1. For each edge, the roots of nodes u and v are in the union-find structure to track connected components. If the root node... and If they differ, they are merged, and the persistence of dimension 0 is recorded as follows. Otherwise, initialize the persistence pairs of dimension 1 as follows: ;

[0051] For dimension 0, the formula for calculating the persistent value is:

[0052] (2);

[0053] Where (b, d) represents a persistent pair, where b is the "birth" time and d is the "death" time. For 0 dimensions, b is usually 0; for 1 dimensions, d is ∞, indicating that the cycle is not closed. persistence0 represents the set of all 0-dimensional persistent pairs. After calculating the persistence value in 0 dimensions, the average persistence value can be obtained, and then normalized to avoid numerical instability. The formula is as follows:

[0054] (3);

[0055] (4);

[0056] in Used to prevent division by zero errors;

[0057] For dimension 1, the number of sustaining pairs needs to be calculated, representing the topological holes in the graph, and the formula is as follows:

[0058] (5);

[0059] If there are no edges in the graph, the default value [0, 0, 0, 0] is returned. These features compactly characterize the persistent cohomological structure of the graph, enhancing the model's ability to incorporate topological signals into gene embeddings.

[0060] Figure 2 This diagram illustrates how persistent graphs are used in the TDA compute framework to preserve the homology information of a given graph through a graph filtering process.

[0061] Figure 3 The diagram illustrates the TDAGENE filtering learning process, which calculates persistent homology from a given graph to obtain TDA feature embeddings. This module effectively extracts global multi-scale information from graph datasets, enhancing the model's ability to extract information from graph data.

[0062] 1.4 Feature Encoder Based on Graph Attention Network

[0063] This invention employs a multi-layered GAT (Gaussian Attention) layer as the feature encoder for TDAGENE. Each GAT layer introduces a self-attention mechanism in the information propagation step, enabling nodes to assign importance weights to their neighbors in the graph. The encoder input is a priori control graph. adjacency matrix and gene expression profile If necessary, normalized node degrees can be selectively concatenated. As an additional feature, the feature dimension concatenated with the additional feature dimension is used as the effective input dimension;

[0064] In each GAT layer, for nodes and his neighbors The attention coefficient is calculated as follows:

[0065] (6);

[0066] in It is a linear transformation matrix. For attention vectors, This indicates a concatenation operation; the negative slope of LeakyReLU is set to... The attention coefficient is normalized by applying softmax to neighboring nodes, and the formula is as follows:

[0067] (7);

[0068] in Let i represent the set of neighbors of node i. The updated embedding is represented as:

[0069] (8);

[0070] This model uses four GAT layers. The outputs of each multi-head attention module are aggregated based on the mean or concatenation. Additionally, each layer introduces residual connections with linear projection to ensure training stability. The final result is a GAT embedding that captures local and contextual gene interactions, as shown in the following formula:

[0071] (9).

[0072] 1.5 Fusion Layer

[0073] This invention employs a gating fusion mechanism to embed GAT. TDA features The ensemble is performed to generate enhanced representations by repeating and adjusting the TDA vector of each node to match the dimension of the GAT embedding, thus forming an expanded topological feature matrix. Subsequently, the GAT embedding is concatenated with the topological features, and a weighted gating vector is generated through a Sigmoid gating unit, as shown in the following formula:

[0074] (10);

[0075] (11);

[0076] in The Sigmoid function is used to adjust the embedding vector using this gate vector, and its formula is as follows:

[0077] (12);

[0078] in The fused output, representing the Hadama product, is further refined by the projection MLP to embed topological information to support downstream branch processing.

[0079] 1.6 Branch-specific multilayer perceptron

[0080] This invention designs a branch-specific multilayer perceptron, and the fused embedding The transcription factor branch and target gene branch are processed independently. Each branch contains three linear layers. Each module sequentially includes linear transformation, layer normalization, LeakyReLU activation, and Dropout regularization. Finally, the final embedding vector of the transcription factor branch is output through a projection layer, as shown in the following formula:

[0081] (13);

[0082] (14);

[0083] The target gene branch is processed using a symmetrical structure, and its calculation process is the same as that of the TF branch:

[0084] (15);

[0085] This dual-branch structure enables specialized feature learning for different regulatory roles of transcription factors and target genes. Finally, LeakyReLU is used as the output layer activation function to ensure the generation of nonlinear and expressive structures, providing support for subsequent prediction of regulatory relationship pairings.

[0086] 1.7 Decoder

[0087] This invention employs an MLP decoder to calculate the final regulatory value, based on transcription factor embedding. and target gene embedding The embeddings of the two branches are concatenated and processed through a three-layer neural network. Each layer includes linear transformation, layer normalization, and LeakyReLU, as shown in the following formula:

[0088] (16);

[0089] Where c represents the integration of TF and the target gene, and I represents the final regulatory value.

[0090] 1.8 Training Optimization

[0091] For training the TDAGENE model, this invention uses binary cross-entropy as the loss function to effectively measure the error between the predicted result and the true label, and utilizes its gradient to efficiently update the network weights, thereby driving continuous optimization of model performance.

[0092] (17);

[0093] In the training of the TDAGENE model, this invention adopts a full-batch gradient descent strategy and introduces reparameterization techniques to achieve effective gradient estimation. The optimization process uses the Adam optimizer with a learning rate of 0.005 and a weight decay coefficient of 0.0005. All key hyperparameters are systematically optimized through grid search.

[0094] Experimental results

[0095] 2.1 Benchmark Dataset

[0096] This invention evaluated the performance of TDAGENE in gene regulatory network (GRN) inference using seven widely used scRNA-seq datasets derived from two human cell lines and five mouse cell lines from BEELINE, including human embryonic stem cells (hESC), human mature hepatocytes (hHEP), mouse embryonic stem cells (mESC), mouse dendritic cells (mDC), and cell lineages from mouse hematopoietic stem cells (mHSC), including erythroid mHSC-E, lymphoid mHSC-L, and granulocyte-monocyte lineage mHSC-GM. Each dataset contained real interaction networks from three different sources: STRING, cell type-specific ChIP-seq, and non-specific ChIP-seq. Detailed information on each dataset and its respective real network is summarized in Table 1.

[0097] Table 1. Statistics of seven single-cell sequencing datasets and training set size for each real network with transcription factors and 500 (1000) most variable genes.

[0098]

[0099] 2.2 Evaluation Strategy

[0100] This invention employs a hold-out validation strategy, treating gene regulatory interactions in real-world networks as positive samples. Positive samples are allocated 8 / 10 to the training set, 1 / 10 to the validation set, and the remainder to the test set. Two different methods are used to select negative samples. For STRING and non-specific ChIP-seq, negative samples are sampled at a 1:1 ratio to ensure a balance between positive and negative samples. For cell type-specific ChIP-seq, a difficult negative sample sampling strategy is used, explicitly defining negative samples as all genes not positively correlated with the given TF. These negative samples are then allocated to the training, validation, and test sets in the same proportion as the positive samples, with the number of negative samples limited by a specified positive-to-negative ratio. For performance evaluation, this invention uses two commonly used prediction metrics: Area Under the ROC Curve (AUROC), representing the balance between true positive and false positive rates at different thresholds; and Precision-Recall (AUPRC).

[0101] 2.3 Inference of gene regulatory network performance in seven benchmark tests

[0102] This invention compares TDAGENE with nine state-of-the-art gene regulatory network inference algorithms, including GRANet, GCLink, GATCL, GNNLink, GENELink, STGRNS, GNE, DeepSEM, and GENIE3. To ensure fairness, this invention uses the same training-test splitting strategy in all methods. This invention ranks TF-gene interactions based on regulatory scores for evaluation.

[0103] To verify the performance advantages of TDAGENE, this invention evaluated TDAGENE and several other methods on datasets derived from seven cell types, including STRING, non-specific, and specific networks. The datasets contained 500 and 1000 highly variable genes. Experimental results showed that the TDAGENE method consistently and stably outperformed existing baseline models such as GRANet and GCLink on the key metric of AUROC, fully validating its superior performance. Figure 4 , Figure 5 and Figure 6 As shown;

[0104] Specifically, compared with Granet, GCLink, GATCL, GNNLink, GENELink, STGRNS, GNE, DeepSEM, and GENIE3, TDAGENE's average performance was improved by 3.08%, 6.31%, 5.45%, 8.22%, 8.44%, 30.57%, 25.61%, 39.18%, and 47.09%, respectively.

[0105] Similarly, Figure 7 , Figure 8 and Figure 9 TDAGENE was evaluated using the AUPRC metric, which demonstrates excellent performance on imbalanced datasets. TDAGENE achieved the best performance on 90% of the datasets, reflecting the model's superior predictive ability. On average, TDAGENE outperformed Granet, GCLink, GATCL, GNNLink, GENELink, STGRNS, GNE, DeepSEM, and GENIE3 by 17.66%, 30.25%, 29.92%, 36.11%, 46.61%, 33.84%, 123.54%, 160.51%, and 170.21%, respectively.

[0106] Figure 4The results show TDAGENE’s scores against nine comparative methods (a total of 10 models) in three real networks for transcription factors (TFs) and AUROC scores for 500 (1000) of the most variable genes. The results are based on a non-specific ChIP-seq dataset and are used to evaluate the performance of each method in identifying transcription factor binding sites and regulatory relationships of highly variable genes.

[0107] Figure 5 The results show TDAGENE's scores against nine other methods in three real-world networks for transcription factors (TFs) and AUROC scores for 500 (1000) of the largest variable genes. The results are based on cell type-specific ChIP-seq datasets and are used to measure the accuracy of each method in predicting gene regulatory networks in specific cell types.

[0108] Figure 6 The results show TDAGENE and nine contrasting methods in three real networks for transcription factors (TFs) and AUROC scores with 500 (1000) maximum variable genes, based on the STRING protein interaction dataset, to evaluate the predictive ability of each method when integrating known interaction information.

[0109] These results further highlight the excellent recognition ability of the TDAGENE model on positive samples. In the inference of single-cell gene regulatory networks, the number of positive examples is much smaller than that of negative examples. The outstanding performance on different datasets consistently demonstrates the reliable potential of TDAGENE in inferring gene regulatory networks.

[0110] Figure 7 The results show the scores of TDAGENE and nine comparative methods (a total of 10 models) on three real networks for transcription factors (TFs) and AUPRC with 500 (1000) maximum variable genes. The results are based on a non-specific ChIP-seq dataset and are used to evaluate the balance between precision and recall of each method in identifying transcription factor binding sites and regulatory relationships of highly variable genes.

[0111] Figure 8 The results show TDAGENE and nine comparative methods in three real networks for transcription factors (TFs) and AUPRC scores with 500 (1000) maximum variable genes. The results are based on cell type-specific ChIP-seq datasets and are used to measure the precision and recall of each method in predicting gene regulatory networks in specific cell types.

[0112] Figure 9The results show TDAGENE and nine comparative methods in three real networks for transcription factors (TFs) and AUPRC scores with 500 (1000) maximum variable genes. The results are based on the STRING protein interaction dataset and are used to evaluate the precision and recall of each method when integrating known interaction information.

[0113] Furthermore, graph-based supervised models, such as GRANet and GCLink, outperform other traditional machine learning and other neural network architectures on two key metrics, demonstrating the growing importance of the graph-based message passing paradigm in the field of gene regulatory networks.

[0114] Compared to traditional machine learning methods, graph neural networks (GRNs) provide a powerful end-to-end learning framework for gene regulatory network inference. Their core advantage lies in their ability to naturally combine the graph structure attributes of biological systems with the powerful representation learning capabilities of deep learning. They can integrate multi-omics data and learn the functional representation of genes in the global network context, thereby more accurately predicting complex and nonlinear regulatory relationships. TDA provides a unique and powerful perspective for gene regulatory network inference. The advantage of TDA is not in accurately predicting a specific regulatory edge, but in revealing the overall structure, stability, and functional modules of GRNs from a global and topological perspective. TDA treats GRNs as a high-dimensional geometric shape, and it is concerned with the global features of this shape, such as connectivity and loop structure, which are not affected by small perturbations. These features often correspond to important biological properties.

[0115] 2.4 Visualization of GRN inference using TDAGENE

[0116] To verify the reliability of TDAGENE in identifying synergistic relationships between genes, this invention visualized gene embeddings using UMAP graphs and performed clustering on the hESC dataset using the Leiden method, such as... Figure 10 As shown, this invention has discovered that the three core transcription factors involved in the pluripotency and self-renewal maintenance of embryonic stem cells—NANOG, SOX2, and POU5F1—aggregate together. A higher number of TDA rings generally indicates a more complex distribution of data points, which can be interpreted as greater heterogeneity in cell states or populations. Figure 11 and Figure 12As shown, from 12h to 24h, the cell population was relatively homogeneous with a less complex topological structure. At 36h, the distribution of biological indicators became more dispersed, and the TDA ring number reached its maximum. This strongly indicates that the cell population was highly heterogeneous at this time, possibly undergoing cell differentiation or state changes, forming a complex and continuous gene interaction expression, which was topologically manifested as a large number of "rings". From 72h to 96h, the dynamic and transitional process may have ended, and the cell population may be converging to a few new and stable states. Therefore, the topological structure is simplified and the ring number decreases. H9, as the control group, may have maintained a relatively simple and stable pluripotent state. In short, the TDA ring number serves as a quantitative indicator of the dynamic changes in intracellular gene expression, and it is closely related to the distribution width and morphology of biological indicators representing cell characteristics.

[0117] Figure 10 The results of UMAP dimensionality reduction visualization of gene embedding vectors output by TDAGENE on the human embryonic stem cell (hESC) dataset are presented to intuitively show the distribution patterns and population structure of genes during dynamic differentiation.

[0118] Figure 11 The dynamic changes of gene regulatory loops (or loop numbers) identified by the TDAGENE method in time series (e.g., from 12h to 96h) are shown, reflecting the temporal evolution of regulatory relationships within cell populations;

[0119] Figure 12 The study showed changes in gene expression levels of three core transcription factors, NANOG, SOX2, and POU5F1, in human embryonic stem cells from 12 to 96 hours, using the H9 cell line as a benchmark to assess the significance of expression changes in the time series.

[0120] Figure 13 A visualization of a real-world network for predicting cell types using TDAGENE is presented. In the diagram, nodes represent the top 30 genes with the highest degree, and edges represent regulatory relationships.

[0121] Figure 14 The degree distribution and main topological properties of the gene regulatory network constructed by TDAGENE on the hESC dataset are shown. The x-axis represents the degree of the network (denoted as k), and the y-axis represents the frequency of the network degree (denoted as P(k)). Both the x-axis and y-axis are logarithmic. The slope of the degree distribution is placed in the upper right corner, and the inset in the lower left corner is the degree distribution of the prior regulatory network.

[0122] Figure 15 It demonstrates the number of transcription factor (TF) regulation events in a real prior regulatory network;

[0123] Figure 16The diagram illustrates the impact of four hyperparameter values ​​on model performance: alpha (a), dropout (b), learning rate (c), and embedding dimension (d). The blue line represents the AUROC value, and the orange line represents the AUPRC value.

[0124] To further evaluate the performance of the TDAGENE inferred gene regulatory network, this invention demonstrates the top 30 genes with a degree greater than 500 inferred GRNs on the cell type-specific dataset hESC 500, such as... Figure 13 As shown in the figure, this invention has found that the inferred gene regulatory networks exhibit scale-free properties, meaning their degree distribution follows a power law, such as... Figure 14 As shown, this is a characteristic of most biological networks. This invention examines the regulatory relationships between eight transcription factors and their ten target genes with the highest regulatory scores. This invention reveals that the reconstructed gene regulatory network can reconstruct the regulatory relationships of key transcription factors such as TFAP2A and JUND. Furthermore, TDAGENE also demonstrates excellent prediction results for transcription factors GTAT4 and KLF5, for which prior data is limited. Figure 15 As shown, JUND plays a regulatory role in genes required for cell cycle progression and is a key regulator of cell cycle processes. Table 2 also shows the top ten gene pairs with the highest prediction scores. These predictions were validated using the publicly available hTFtarget database. The interactions between seven gene pairs have been confirmed, while the remaining three pairs, although not currently verifiable, may represent real gene interactions. Comprehensive proteomics, phosphoproteomics, and generalization analyses of tumors and adjacent normal liver tissues from liver cancer patients showed high expression of COL4A1 in patients with poor disease-free survival. Studies have shown that the maternal transcription factor OTX2 regulates embryonic genome activation and early embryogenesis through epigenetic mechanisms. These findings demonstrate that TDAGENE not only has excellent predictive performance for confirmed gene pairs but also holds great potential for predicting new gene pairs, opening up entirely new perspectives for researchers to discover gene regulatory relationships.

[0125] Table 2. Top 10 gene pairs predicted by TDAGENE on the hESC dataset

[0126]

[0127] 2.5 Experimental Parameter Settings

[0128] In gene regulatory network inference, the four hyperparameters of the TDAGENE model play a crucial role: learning rate, embedding dimension, dropout rate, and alpha in the GAT layer. To determine the optimal configuration, this invention employs a grid search strategy, adjusting each hyperparameter individually while keeping other parameters constant. Hyperparameter analysis in this invention uses the hHSC dataset from the STRING standard network, and performance is measured using average AUROC and AUPRC across 10 experiments. Detailed results for each hyperparameter pair are as follows: Figure 16 As shown;

[0129] The parameter alpha adjusts the LeakyReLU activation function, affecting the slope of the negative half-axis. This invention tested the alpha value within the range of 0 to 0.9, setting a step size of 0.05. Figure 16 As shown in Figure (a), the results indicate that TDAGENE achieves the best performance on both important performance metrics when the alpha value is set to 0.4. The choice of dropout is a crucial trade-off. Using dropout during training to randomly discard neurons can lead to the loss of key information and prevent the network from learning effectively. Conversely, a dropout rate that is too low results in very weak regularization and cannot effectively prevent overfitting. This invention tested dropout rates from 0 to 0.9 with a step size of 0.05. Experiments show that the model performs best when the dropout rate is 0.1. Figure 16 As shown in Figure (b), the learning rate is a key hyperparameter affecting the convergence speed and quality of deep learning models. In the experiments of this invention, the learning rate is set to 5e. -4 At that time, TDAGENE achieved the best performance on AUROC and AUPRC, such as Figure 16 As shown in Figure (c), this invention tested the gene embedding dimensions of 8, 16, 32 and 64. The model performed best at an embedding dimension of 32. Higher dimensions lead to more complex and redundant calculations, which do not show better performance but increase computational costs. Lower dimensions do not fully express genes, and the lower performance indicates that low-dimensional gene expression is not conducive to the inference of gene regulatory networks.

[0130] The hyperparameter settings of this invention fine-tuned the TDAGENE model, significantly improving its performance in gene regulatory network inference tasks through fine-tuning. These findings highlight the key role of parameter optimization in achieving advanced performance in complex neural network models.

[0131] 2.6 Ablation Study

[0132] To verify whether the submodules of TDAGENE contribute to the inference performance of gene regulatory networks and to quantify its improvement, this invention conducted a series of ablation experiments by progressively removing or replacing key components of the model and analyzing their impact on performance. Specifically, this invention constructed the following four variants, as shown in Table 3:

[0133] 1) TDAGENE without TDA: Completely removes TDA computation and fusion mechanisms, and only uses GAT embedding to directly input MLP branches;

[0134] 2) TDAGENE without fusion: TDA features and GAT embeddings are directly concatenated without using a fusion mechanism;

[0135] 3) TDAGENE with no residual joins: Retains four layers of GAT but removes residual joins;

[0136] 4) TDAGENE has no independent branch: TF and target genes are merged into a shared MLP branch, and a unified embedding is output before decoding;

[0137] Each variant used the same hyperparameter settings as TDAGENE, and was tested ten times on a dataset of seven cell types, with mean AUROC and AUPRC used as metrics for evaluation. Figure 17 and Figure 18 As shown, removing any key module leads to a significant decrease in model performance. Notably, removing the TDA module causes the most significant performance drop, with AUROC decreasing by 3.71% and AUPRC by 32.87%. TDA systematically identifies and analyzes stable and synergistic regulatory relationships between genes. In the topology, 1-dimensional "holes" typically correspond to feedback loops, which are the core of the regulatory network and determine the system's stability and oscillatory behavior, such as biological clocks and cell fate determination. TDA can automatically and globally scan the entire network to identify all these potential feedback structures. The results show that this is effective for identifying positive samples in imbalanced datasets. Removing residual connections leads to a 16.33% decrease in AUROC and an 18.99% decrease in AUPRC, while removing TF-target gene branching leads to a 5.98% decrease in AUROC and a 20.75% decrease in AUPRC. In gene regulatory network inference, branch specialization helps to better identify and establish the relationships between genes. These results indicate that the TDA module, fusion module, residual connections in GAT, and TF-target gene branching all significantly contribute to the performance improvement of TDAGENE.

[0138] Table 3. Variations and modifications of TDAGENE

[0139]

[0140] A checkmark indicates that the module is in the model, while a cross indicates that it is not.

[0141] 2.7 Effects of Filter Scale Variation

[0142] Figure 17 The ablation experiments of the four main components of TDAGENE, four variant models, and the AUROC values ​​of TDAGENE on seven different cell type datasets are presented. The average performance degradation of each variant model is 3.71%, 10.59%, 16.33%, and 5.98%, respectively.

[0143] Figure 18 The ablation experiments of the four main components of TDAGENE, four variant models, and the AUPRC values ​​of TDAGENE on seven different cell type datasets are presented. The average performance degradation of each variant model was 32.87%, 14.28%, 18.99%, and 20.75%, respectively.

[0144] A significant advantage of TDAGENE lies in its embedding of TDA features into gene expression embeddings. TDA provides global structural insights by transforming a complex gene interaction network into a persistent graph representation. The choice of filter scale directly affects the invention's assessment of persistence. At too small a filter scale, only very close points are connected, resulting in a large number of isolated small connected components, making it impossible to see the overall structure. As the scale increases, some connected components merge, and some loop structures may form. When the scale is large enough, all points eventually connect into a large connected component, filling all holes. However, excessively large scales lead to very complex simple complexes, with an explosive increase in the number of vertices, edges, and faces, resulting in a sharp rise in computational costs. This invention tested filter scales from 0 to 2 with a step size of 0.5, as shown in Table 4. The model achieved the best performance when the filter scale was set to 1.0. Excessively high filter scales increased computational costs without providing a significant improvement.

[0145] Table 4. Effect of scaling factor on filter scale

[0146]

[0147] Single-cell RNA sequencing technology provides an unprecedented perspective for analyzing the heterogeneity of gene regulation at the cellular level. However, accurately inferring gene regulatory networks from high-dimensional, sparse, and noisy scRNA-seq data, especially capturing their intrinsic global structure and complex nonlinear relationships, remains a core challenge for computational biology. Traditional unsupervised methods are often accompanied by high false positives, while existing supervised models, although improving accuracy, are mostly limited to local node interactions and are difficult to understand the stability and functional modules of the regulatory system from the perspective of the overall topology.

[0148] To address this, this invention proposes TDAGENE, an innovative framework that integrates topological data analysis with graph attention networks. The core concept of this method is to treat the gene regulatory network as a geometric shape with rich topological features. By using TDA techniques such as persistent cohomology, its global multi-scale structure, such as connected components and feedback loops, is quantified, thereby compensating for the shortcomings of traditional graph learning methods in capturing global features. TDAGENE adaptively combines these topological insights with the local gene embeddings learned by GAT through a gated fusion mechanism, and uses a branch-specific multilayer perceptron to perform specialized characterization learning of transcription factors and target genes, ultimately achieving accurate and interpretable inference of gene regulatory relationships.

[0149] Systematic evaluations on multiple benchmark datasets demonstrate that TDAGENE significantly and consistently outperforms a range of state-of-the-art methods on key metrics such as AUROC and AUPRC, fully validating its superior performance. Importantly, in cell-type-specific networks with high imbalance, TDAGENE achieves particularly outstanding improvements on AUPRC due to its ability to capture global topology. Ablation studies further quantify the contributions of each core component, confirming that TDA features, gated fusion, residual connections, and bi-branch architecture are all crucial to the model's success. Case studies show that TDAGENE can not only cluster known pluripotency factors and identify topological dynamics during cell state transitions, but also predict novel regulatory relationships validated by databases, demonstrating its enormous potential for discovering new biological knowledge.

[0150] This invention provides a novel and powerful paradigm for inferring single-cell gene regulatory networks by injecting a global perspective of topological data analysis into graph representation learning. It not only significantly improves prediction accuracy, but also provides new theoretical and tool support for understanding cell fate determination, heterogeneity, and disease mechanisms at the system level.

[0151] The above embodiments are illustrative of the present invention and are not intended to limit the present invention. Any simple modifications to the present invention are within the scope of protection of the present invention.

Claims

1. A gene regulation inference method guided by topological data analysis for gene network embedding, characterized in that: Includes the following steps: S1. Topological Data Analysis Feature Calculation: Calculate TDA features from gene expression data and prior knowledge to capture the global topological characteristics of the gene regulatory graph. The persistence value of the Betti number reflects the change of topological features during the graph filtering process. The persistence value is obtained by using a filtering function. The persistence process of the graph is constructed using persistent graphs and persistent barcodes. The TDA feature embedding is obtained using persistent barcodes. S2. Feature Encoder Based on Graph Attention Network: A multi-layer GAT layer is used as the feature encoder of TDAGENE. Each GAT layer introduces a self-attention mechanism in the information propagation step, and assigns importance weights to neighbors in the graph through nodes. There are four GAT layers. The output of each multi-head attention module is aggregated according to the mean or splicing. Each layer introduces residual connections with linear projection to generate GAT embeddings that capture local and contextual gene interaction relationships. S3, Fusion Layer: A gated fusion mechanism is used to integrate GAT embeddings with TDA features to generate enhanced representations. The TDA vector of each node is repeated and adjusted to match the dimension of the GAT embedding, forming an expanded topological feature matrix. Then, the GAT embeddings and topological features are concatenated, and a weighted gate vector is generated through a Sigmoid gate unit. The embedding vector is adjusted through this gate vector. The fusion output is further refined by a projection MLP so that the embedding contains topological information. S4, Branch-Specific Multilayer Perceptron: The fused embeddings are processed through independent transcription factor branches and target gene branches; each branch contains three linear layers, and each module sequentially includes linear transformation, layer normalization, LeakyReLU activation, and Dropout regularization. The final embedding vector of the transcription factor branch is output through a projection layer, and the target gene branch is processed using a symmetrical structure. Specific features are learned for different regulatory roles of transcription factors and target genes, and LeakyReLU is used as the activation function of the output layer. S5, Decoder: The final regulatory value is calculated using an MLP decoder. Based on the transcription factor embedding and the target gene embedding, the embeddings of the two branches are spliced ​​together and processed through a three-layer neural network. Each layer contains linear transformation, layer normalization and LeakyReLU to decode and predict the regulatory links between gene pairs. S6. Training Optimization: Use binary cross-entropy as the loss function to measure the error between the predicted result and the true label, and use its gradient to update the network weights to drive the continuous optimization of model performance. We employ a full-batch gradient descent strategy, introduce reparameterization techniques for gradient estimation, use the Adam optimizer, set the learning rate to 0.005 and the weight decay coefficient to 0.0005, and optimize all key hyperparameters through grid search.

2. The method of gene regulation inference guided by topological data analysis for gene network embedding of claim 1, wherein: In step S1, a graph structure is constructed using prior knowledge, and a filtering function is used to identify global topological features of the graph from the graph structure, including connected components and one-dimensional holes. Persistent values ​​reflect changes in global topological features, and persistent value graphs and persistent value barcodes are used to obtain the embedding of global topological features.

3. The method of gene regulation inference guided by topological data analysis-driven gene network embedding of claim 1, wherein: The step S1 extracts and represents the topological features of the graph through the graph filtering process, models by using a neural network, trains the filtering function for the graph learning task to optimize the filtering process, the filtering function is based on nodes and edges, obtains the TDA feature embedding by calculating persistent homology of a given graph, and extracts global multi-scale information of the graph dataset.

4. The method of gene regulation inference guided by topological data analysis for gene network embedding of claim 1, wherein: In the step S2, the feature encoder takes the sparse adjacency matrix representing the prior regulatory relationship and the gene expression matrix as input, splices the normalized node degree as additional features, and takes the feature dimension of the spliced input dimension and additional feature dimension as the effective input dimension.

Images