Gene regulation and control network construction method, system, equipment and medium

Abstract

The invention belongs to the related technical field of gene regulation and control networks, and provides a gene regulation and control network construction method, system, equipment and medium in order to overcome the defects of the existing gene regulation and control network in the aspects of interpretability, authenticity and the like. Based on a graph neural network, efficient gene embedding representation learning is realized to capture a global dependency relationship; and then based on different types of regulation and control relationships, analyzing a linear causal relationship and a nonlinear causal relationship, integrating the two channels to obtain a confidence coefficient matrix of the gene, and further constructing a gene regulation and control network. The gene regulatory network constructed by the invention can reveal the regulatory relationship between genes, the interpretability, authenticity and credibility of the gene regulatory network are improved, and a more accurate and reliable tool is provided for analyzing a complex biological system subsequently.

Description

Technical Field

[0001] The present invention belongs to the technical field related to gene regulatory networks, and in particular relates to a gene regulatory network construction method, system, equipment and medium. Background Art

[0002] The statements in this section merely provide background information related to the present invention and do not necessarily constitute prior art.

[0003] As the core mechanism for gene expression regulation in biological systems, the gene regulatory network (GRN) is of great significance for understanding cell development, differentiation, physiological function, and the occurrence and development of diseases. From the perspective of feedback regulation mechanism, the gene regulatory network is a complex and dynamic system composed of interactions between chromatin, transcription factors (TFs), and genes. These interactions form a complex regulatory circuit that determines gene expression patterns and cell functions. Each cell type or cell stable state is composed of a specific set of transcription factors and their target genes (TGs).

[0004] In recent years, with the rapid development of high-throughput sequencing technology, a large number of experimental techniques have been applied to the study of gene regulatory networks. For example, the emergence of single-cell sequencing technology has provided an effective method to study the heterogeneity of gene regulatory networks. Single-cell sequencing can collect omics data from individual cells, providing insights into the stochastic nature / probability of gene expression and the underlying regulatory programs. In addition, with the development of computational biology, an increasing number of computational methods are being used to infer gene regulatory networks from large-scale biological data. These methods not only help to better understand the structure and function of gene regulatory networks but also provide new ideas and methods for disease diagnosis and treatment.

[0005] However, existing gene regulatory network inference methods are mostly based on correlation analysis. Gene regulatory networks are essentially a collection of causal relationships between transcription factors and target genes. This results in regulatory networks constructed based on correlations lacking in interpretability, authenticity, and credibility, and they are unable to accurately reveal the regulatory relationships between genes. Furthermore, when faced with real biological data, the high complexity and massive scale of gene regulatory networks pose significant challenges to traditional causal structure learning methods that rely solely on information between pairs of nodes.

[0006] Therefore, how to rely on real observational data to construct an accurate and reliable causal gene regulatory network remains an urgent problem to be solved. Summary of the Invention

[0007] To overcome the deficiencies of the above-mentioned prior art, the present invention provides a method, system, device and medium for constructing a gene regulatory network. The constructed gene regulatory network can reveal the regulatory relationship between genes, thereby improving the interpretability, authenticity and credibility of the gene regulatory network.

[0008] In order to achieve the above object, the present invention adopts the following technical solutions: In a first aspect, the present invention provides a method for constructing a gene regulatory network, comprising: According to the gene expression matrix, the linear and nonlinear relationships between genes were determined using the Pearson correlation coefficient method and mutual information, and the first and second gene co-expression networks were constructed respectively; Based on the first gene co-expression network and the gene expression matrix, and the second gene co-expression network and the gene expression matrix, respectively, a graph convolutional network model is used to capture global dependencies for gene embedding representation learning to obtain a first gene embedding matrix and a second gene embedding matrix; Based on the first gene embedding matrix and the second gene embedding matrix, according to different types of regulatory relationships, a causal relationship score is calculated using an additive noise model to obtain a first confidence matrix and a second confidence matrix, and the first confidence matrix and the second confidence matrix are combined to obtain a causal confidence matrix of the gene; A gene regulatory network is constructed based on the causal confidence matrix of the gene and the gene expression matrix.

[0009] In a second aspect, the present invention provides a gene regulatory network construction system, comprising: A construction module is configured to: determine the linear relationship and the nonlinear relationship between genes according to the gene expression matrix using the Pearson correlation coefficient method and mutual information, and respectively construct a first gene co-expression network and a second gene co-expression network; an embedding representation module configured to: based on the first gene co-expression network and the gene expression matrix, and the second gene co-expression network and the gene expression matrix, respectively use a graph convolutional network model to capture global dependencies and perform gene embedding representation learning to obtain a first gene embedding matrix and a second gene embedding matrix; An embedding representation module is configured to: calculate a causal relationship score based on the first gene embedding matrix and the second gene embedding matrix according to different types of regulatory relationships using an additive noise model to obtain a first confidence matrix and a second confidence matrix, and synthesize the first confidence matrix and the second confidence matrix to obtain a causal confidence matrix of the gene; The gene regulatory network module is configured to construct a gene regulatory network based on the causal confidence matrix of the gene and the gene expression matrix.

[0010] In a third aspect, the present invention provides an electronic device comprising a memory and a processor, and computer instructions stored in the memory and executed on the processor, wherein the computer instructions, when executed by the processor, perform the method described in the first aspect.

[0011] In a fourth aspect, the present invention provides a computer-readable storage medium for storing computer instructions, wherein when the computer instructions are executed by a processor, the method described in the first aspect is performed.

[0012] One or more of the above technical solutions have the following beneficial effects: In the present invention, gene co-expression networks are constructed respectively by using the Pearson correlation coefficient and mutual information, and efficient gene embedding representation learning is achieved based on graph neural networks to capture global dependencies; then, based on different types of regulatory relationships, causal relationship scores are calculated using an additive noise model, linear causal relationships and nonlinear causal relationships are analyzed, and the confidence matrices obtained by linear and nonlinearity are integrated to obtain the causal confidence matrix of the gene, thereby constructing a gene regulatory network. The present invention provides a new paradigm for the inference of causal gene regulatory networks. The constructed gene regulatory network can reveal the regulatory relationships between genes, improve the interpretability, authenticity and credibility of the gene regulatory network, and subsequently provide a more accurate and reliable tool for analyzing complex biological systems.

[0013] Advantages of additional aspects of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] The accompanying drawings, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0015] Figure 1 Schematic diagram of the existing causal learning GNN framework; Figure 2 Schematic diagram of the pairwise causal discovery framework combined with clustering graph neural network proposed in Example 1 of the present invention; Figure 3 Schematic diagram of the gene embedding representation learning implementation process in Example 1 of the present invention; FIG4 (a) is a schematic diagram of the gene regulatory network prediction performance of different methods on the hESC-Specific-500 dataset and the hESC-CellType-1000 dataset under different evaluation indicators in Example 1 of the present invention; FIG4( b ) is a schematic diagram of the gene regulatory network prediction performance of different methods on the mDC-CellType-500 dataset and the mDC-CellType-1000 dataset under different evaluation indicators in Example 1 of the present invention; FIG5 (a) is a schematic diagram showing the comparison results of gene proportions in knockout data in the TCGAS gene library in Example 1 of the present invention; FIG5( b ) is a schematic diagram showing the comparison results of gene proportions in knockout data in the Rembrandt gene library in Example 1 of the present invention; FIG5 (c) is a schematic diagram showing the comparison results of the gene proportions in the knockout data of the CGGA_301 gene library in Example 1 of the present invention; FIG6 (a) is a schematic diagram showing the effect of TCGA_GRHL3 expression on patient survival in Example 1 of the present invention; FIG6 (b) is a schematic diagram showing the effect of Rembrandt_GRHL3 expression on patient survival in Example 1 of the present invention; FIG6( c ) is a schematic diagram showing the effect of CGGA_301_GRHL3 expression on patient survival in Example 1 of the present invention. DETAILED DESCRIPTION

[0016] It should be noted that the following detailed descriptions are exemplary and intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0017] It should be noted that the terms used herein are for describing particular embodiments only and are not intended to limit the exemplary embodiments according to the present invention.

[0018] In the absence of conflict, the embodiments of the present invention and the features thereof may be combined with each other.

[0019] Gene regulatory network inference is a key research direction for unraveling the complex regulatory mechanisms of biological systems. Its core lies in accurately identifying the regulatory relationships between transcription factors and target genes from gene expression data. Current research methods in this field are diversified and can be broadly categorized into two main approaches: correlation-based and causality-based. The former relies on statistical association analysis to construct regulatory networks between genes, forming a rich cluster of methods within a data-driven modeling framework. The latter aims to uncover the causal nature of regulatory relationships, thereby overcoming the inherent limitations of correlation analysis.

[0020] Currently, correlation-based GRN inference methods mainly fall into the following categories: information-theoretic, regression-based, probabilistic graphical model-based, and machine learning and deep learning-based methods. Information-theoretic methods infer gene regulatory relationships by calculating statistical correlations in gene expression data. For example, the CMI2NI method uses mutual information to quantify gene dependencies. However, this method has difficulty distinguishing between direct and indirect regulation and is insufficiently sensitive to nonlinear relationships between genes. Therefore, it cannot accurately reflect complex biological regulatory mechanisms and has significant limitations. Regression-based methods typically regress the gene expression of a target gene onto a set of transcription factors and use the regression coefficient as the regulatory strength between the transcription factor and the target gene. GRNBoost is a typical example of this method. This method has some effectiveness in handling linear relationships, but when faced with complex real-world regulatory networks, it cannot capture all regulatory relationships, resulting in incomplete constructed gene regulatory networks and selection bias. Methods based on probabilistic graphical models (PGMs), such as Bayesian networks, represent intergene dependencies. BoolODE is a typical example of this approach, which transforms Boolean models into ordinary differential equations (ODEs) to simulate the dynamic behavior of gene regulatory networks. While theoretically capable of handling complex network structures, these methods suffer from complex parameter estimation, high computational complexity, and high computational cost, making them difficult to handle intricate and large real-world regulatory networks. Machine learning and deep learning-based methods, such as the GENIE3 method, use regularization techniques to infer gene regulatory networks. However, this method requires a large amount of labeled data, resulting in poor prediction accuracy for unlabeled data and high data requirements. The ambiguity of the "group decision-making" of ensemble tree models also leads to poor model interpretability. The DeepGRN method, on the other hand, utilizes deep neural networks to infer gene regulatory networks. It exhibits advantages in handling large datasets, overcoming the difficulties of traditional methods in handling multivariate data. However, this method requires a large amount of training data and is prone to overfitting on small datasets. Furthermore, the "black-box" nature of deep neural networks also contributes to poor model interpretability.

[0021] Because GRNs are essentially networks composed of numerous causal relationships from TFs to TGs, GRN inference methods based on causality are theoretically more realistic. From the perspective of data type, causal methods can be divided into two categories: those based on time-tagged / series data and those based on time-untagged data. Causal discovery methods based on time series data, such as Granger causality (GC), use lagged values ​​of a time series to predict future values, thereby determining causal relationships between variables. If the past value of one series significantly improves the prediction accuracy of another, Granger causality is considered to exist between the two. However, this method relies heavily on the prior conditions of the time series data, namely, the data must have clear time labels, that is, there must be a temporal relationship between the series. However, collecting such data typically requires elaborate and complex experiments and continuous monitoring, which are costly and technically limited. This makes it extremely difficult to obtain biometric data with precise time labels, which greatly restricts the application of GC methods. Causal discovery methods based on time-untagged data can be divided into constraint-based methods and functional causal models. The PC algorithm is a classic constraint-based causal discovery method. It first constructs a causal skeleton and then determines the causal direction using conditional independence constraints. Numerous improved causal structure learning methods have been derived from the PC algorithm, such as DCI, IDA, and fast causal inference algorithms. While these methods have made some progress in handling large-scale variables, their computational efficiency and high cost remain unresolved. Functional causal models, such as NOTEARS, combine structural equation models (SEM) to learn directed acyclic graphs (DAGs) from observational data. However, the DAG assumption requires the absence of any feedback loops in the network, which is unrealistic in gene regulatory networks. Regulatory relationships in biological systems are often complex and mutually influential. This assumption is too strong for GRN reasoning and cannot be directly applied to GRN inference tasks.

[0022] GRINCD is an innovative gene regulatory network inference method that combines graph representation learning with causal asymmetry learning to predict regulatory relationships between genes using ANM. This approach considers both linear and nonlinear factors, tackling complex feedback loops in GRNs and the inapplicability of strong DAG assumptions, thereby improving inference accuracy. Leveraging the excellent learning capabilities of graph neural networks for graph-structured data, the GraphSAGE model is used to sample neighboring nodes and aggregate information to learn node embeddings that reflect the characteristics of neighboring nodes, thereby enriching the information of each node. Furthermore, ANM is used to compute pairwise causal regulatory relationships for all TF-TG and TF-TF pairs based on node embeddings, relaxing the DAG assumption and making the GRN construction more reasonable. However, the GraphSAGE model suffers from low memory efficiency when processing large-scale graph data and also suffers from computational efficiency limitations, particularly during the training of deep networks. GraphSAGE's neighbor sampling method incurs significant computational overhead and slows convergence when used with large-scale, high-dimensional, and complex real-world biological observational data. In addition, GraphSAGE suffers from the vanishing gradient problem when training deep networks, which makes it difficult to effectively optimize the underlying network parameters close to the input layer. The model cannot fully extract basic features, thereby reducing the overall performance of the model.

[0023] Unlike GraphSAGE, ClusterGCN is a novel graph convolutional network (GCN) algorithm. It utilizes graph clustering structures for stochastic gradient descent (SGD) training, a feature that significantly improves memory utilization and computational efficiency. During training, ClusterGCN samples a set of nodes associated with a dense subgraph at each step and restricts the neighborhood search to that subgraph, effectively reducing the number of node embeddings that need to be stored. This design significantly reduces memory usage and greatly improves computational efficiency when processing large-scale graph data. It also enables the training of deeper networks without excessive time and memory overhead, improving overall model performance and, in turn, prediction accuracy.

[0024] Example 1 Based on the above related work, this embodiment discloses a method for constructing a gene regulatory network, including: Step 1: Based on the gene expression matrix, the linear and nonlinear relationships between genes are determined using the Pearson correlation coefficient method and mutual information, and the first gene co-expression network and the second gene co-expression network are constructed respectively; Step 2: Based on the first gene co-expression network and gene expression matrix, as well as the second gene co-expression network and gene expression matrix, a graph convolutional network model is used to capture global dependencies and perform gene embedding representation learning to obtain the first gene embedding matrix and the second gene embedding matrix; Step 3: Based on the first gene embedding matrix and the second gene embedding matrix, according to different types of regulatory relationships, the causal relationship score is calculated using the additive noise model to obtain the first confidence matrix and the second confidence matrix, and the first confidence matrix and the second confidence matrix are combined to obtain the causal confidence matrix of the gene; Step 4: Construct a gene regulatory network based on the gene causal confidence matrix and gene expression matrix.

[0025] This example proposes a framework that combines a graph neural network (GNN) and a causal discovery model, aiming to construct a gene regulatory network with causal relationships based on the expression data of transcription factors and their target genes. Figure 1 This example demonstrates a general causal learning (GNN) framework that first performs causal relationship learning and then uses GNNs for information aggregation. This example first leverages the powerful graph-structured data processing capabilities of graph neural networks to extract embedding information for each node from the underlying gene co-expression network, maximizing the ability to capture inter-node dependencies. Furthermore, a method based on the additive noise model (ANM) is proposed to compute pairwise causal regulatory relationships for all TFs-TGs and TFs-TFs. Taking into account both linear and nonlinear correlations in the data, different adaptive methods are combined to construct linear and nonlinear processes, ultimately integrating them to form a complete causal gene regulatory network. Subsequent validation is conducted on two public datasets, and comparisons with 13 widely used gene regulatory network inference methods demonstrate the superior performance of the proposed framework. Furthermore, a combination of dry and wet experiments, using real glioblastoma (GBM) clinical data, test the proposed framework's ability to construct GRNs using complex and high-dimensional real-world biological data.

[0026] The following is a detailed description of a gene regulatory network construction method proposed in this embodiment: First, this embodiment defines the gene regulatory network (GRN) inference task as follows: GRN is considered as a weighted directed bipartite graph ,in, Transcription factors gather, Target gene gather; yes There are directed edge sets in , including all possible and Relationship pairs, especially note that since the edges between target genes do not belong to the discussion scope of GRN, Relationship does not belong to ; Is a causal confidence matrix with a size of ,in, Corresponding to Time Edge or Time Edge This embodiment assumes that there is no self-regulation in GRN, that is, hour In addition, Indicates a size of The gene expression matrix of behavior samples, listed as genes, among which This embodiment will be based on 、 as well as As input, we finally get .

[0027] In step 1 of this example, an initial gene co-expression network is first constructed to provide prior information to guide gene representation learning. Edges in this network represent correlations between genes. The implementation process is shown in Algorithm 1. In the complex regulatory networks of biological systems, gene expression data exhibit a multivariate correlation pattern that interweaves linear and nonlinear correlations. Linear correlations encompass both positive and negative correlations, such as the coordinated expression of the proto-oncogene MYC and the cell cycle protein Cyclin D1, and negative correlations, such as the negative feedback regulation between the tumor suppressor gene TP53 and MDM2. Nonlinear correlations encompass various forms, including quadratic functions, exponentials, logarithms, and periodicity, corresponding to regulatory mechanisms such as threshold effects, signal amplification, saturation responses, and biological rhythms.

[0028] When constructing a gene co-expression network, this embodiment combines the Pearson Correlation Coefficien (PCC) and mutual information (MI) to construct a first gene co-expression network and a second gene co-expression network based on the linear relationship and nonlinear relationship between genes, respectively. Among them, the Pearson correlation coefficient method accurately quantifies the strength and direction of the linear relationship by calculating the ratio of covariance to standard deviation, but lacks the ability to process nonlinear dependence; mutual information is based on the principles of information theory and can capture various linear and nonlinear dependencies, but lacks directional information. The combination of the two can achieve complementary advantages. Mutual information broadens the coverage of association patterns, and the Pearson correlation coefficient method supplements the linear regulation direction, which is more in line with the real data generation process of biological systems and provides a comprehensive and accurate mathematical description for the construction of gene co-expression networks.

[0029] In this embodiment, let is the correlation coefficient matrix, is the correlation coefficient matrix constructed based on the Pearson correlation coefficient method , or the correlation coefficient matrix constructed based on mutual information ,make Set a threshold for something The adjacency matrix of the gene co-expression network generated by To pass Determine the constructed adjacency matrix , or through Determine the constructed adjacency matrix .

[0030] In this embodiment, select As the basic threshold, The value of is shown in formula (1): (1)

[0031] It should be noted that a gene co-expression network typically refers to a graph network consisting of points and lines. The adjacency matrix of a gene co-expression network is a two-dimensional matrix vector, where the row and column names represent points in the graph, and the values ​​indicate whether an edge exists. The adjacency matrix of a co-expression network is typically used in calculations, and the term co-expression network is often used in presentations.

[0032] Traditional causal discovery methods usually use gene expression matrices directly to identify causal regulatory relationships between genes in a system. In this case, the effects of other possible regulatory factors on the two genes under study are often ignored. Figure 3As shown, in order to solve this problem, in step 2 of this embodiment, the adjacency matrix of the gene co-expression network is used And the normalized gene expression matrix As input, a three-layer Cluster-GCN model is used to summarize the effective feature information of each gene and its adjacent genes in the gene co-expression network, and a binary classifier is used for verification.

[0033] Existing work uses the GraphSAGE model to obtain the embedded representation of nodes. It adopts a random walk strategy to select positive and negative links from the gene co-expression network as input, and uniformly samples a set of fixed-size neighborhoods through the GraphSAGE model. In this case, the embedding calculation of each loss term is , where the number of neighbors sampled by each node in each layer is expressed as , the number of model layers is expressed as , that is, when processing large sample sets, the time complexity increases exponentially with the number of model layers according to the sampling size, but the full neighborhood set is not used in this process, which will also reduce the accuracy of gradient estimation.

[0034] This example uses the Cluster-GCN model to obtain node embedding representation. First, for the gene co-expression network , dividing its gene nodes into Group, that is ,in, By The nodes in the partition.

[0035] From this, we can get Sub-networks: (2) Among them, each edge Only node sets The links between them.

[0036] When reorganizing the nodes, the adjacency matrix of the gene co-expression network will be divided into sub-matrices.

[0037] (3) (4) Each diagonal block is a size of The adjacency matrix contains the adjacency matrix from the co-expression network inner edge; It is a subnetwork The adjacency matrix of Contains two partitions and The edge between is the adjacency matrix of the gene co-expression network The matrix consisting of all off-diagonal blocks of .

[0038] Similarly, we can use the The normalized gene expression matrix Divide into multiple sub-matrices , and the label matrix Divided into ,in, and Respectively by The node features and labels in .

[0039] Using the above method, the objective function of the Cluster-GCN model is decomposed into different clusters. Represents the normalized , the final gene embedding matrix becomes: (5) Among them, the activation function In this embodiment, it is set as the ReLU function, is a learnable parameter that is optimized during training by back propagation and represents the The feature transformation matrix of the layer.

[0040] because With a block diagonal form, the loss function It can also be decomposed into Corresponding to the label matrix The Elements: (6) At this point, Cluster-GCN is based on the decomposition form in formula (5) and formula (6). In each step, a cluster is first Sampling is performed, and then according to The gradient of SGD is updated, which only requires the subgraph of the current batch 、 and as well as ,in, This is used to represent the embedding matrix calculation model shown in Formula (5), namely the Cluster-GCN model. Only one block of the matrix product forward and backward propagation in Formula (5) is required during execution, which is easier to implement than the neighborhood search process in the GraphSAGE model. The specific implementation process is shown in Algorithm 2.

[0041]

[0042] In Cluster-GCN, a graph clustering algorithm is used to partition the network. Specifically, the Metis clustering algorithm is used to partition the gene co-expression network. This method can minimize the number of inter-cluster links and thus maximize the embedding utilization.

[0043] Compared to graph neural network methods such as GraphSAGE and VR-GCN, Cluster-GCN has a more linear rather than exponential time complexity and a smaller space complexity. When faced with complex and large inputs such as bioinformatics data, it can complete training tasks with less time consumption and lower memory requirements.

[0044] The detailed time and space complexity of different GCN training algorithms are summarized below: Time complexity of the GCN model: ; Space complexity of the GCN model: ; Time complexity of the Vanilla SGD model: ; Space complexity of the Vanilla SGD model: ; Time complexity of the GraphSAGE model: ; Space complexity of the GraphSAGE model: ; Time complexity of the FastGCN model: ; Space complexity of FastGCN: ; Time complexity of VR-GCN model: ; Space complexity of VR-GCN model: ; Time complexity of the Cluster-GCN model: ; The space complexity of the Cluster-GCN model: ; in, represents the feature vector dimension, which in this embodiment represents the number of samples; represents the average degree of each node, that is, the average number of edges connected to each node; Indicates the batch size; Represents the adjacency matrix The number of non-zero elements of is the total number of edges in the graph, L represents the number of model layers, represents the number of neighbors sampled at each layer of each node, V represents the set of gene nodes, and |V| represents the number of nodes.

[0045] After completing the graph representation learning, this embodiment can obtain a new gene embedding matrix That is, the first gene embedding matrix or the second gene embedding matrix. In , the neighbor node information related to the current gene is aggregated into itself, enriching its own node features.

[0046] In step 3 of this embodiment, the causal relationship is analyzed based on the additive noise model, which is a model that uses nonlinear relationships and Gaussian noise to identify causal directions and is based on the gene embedding matrix. , calculate the causal tendency score between two genes to judge the strength of the causal relationship. Note that GRN in the strict sense (only transcriptional regulation) only exists as well as There are no two types of regulatory relationships. regulatory relationships, as such relationships are usually not directly achievable through transcription.

[0047] ANM assumes that there is a nonlinear relationship between causal variables and the existence of independent additive Gaussian noise. The nonlinear relationship conforms to the complex laws of GRN. Introducing Gaussian noise is also a common method to simulate complex disturbances. For the two different types of regulatory relationships in GRN, the traditional network construction method calculates as well as The propensity scores of all bidirectional relationships are sorted to determine the strength of the relationship. This calculation method will generate a lot of unnecessary calculations when processing the GRN with only transcriptional regulation mentioned in this embodiment, and the calculation cost is high. In this embodiment, different calculation processes are designed for different types of regulatory relationships: Relationship pairs, respectively, two different directions of regression and respectively get the regression residuals to achieve ANM, then by verifying the independence between the residuals and the independent variables to calculate the causal tendency scores in the two directions, and for If the relationship is correct, only calculation in one direction is required.

[0048] According to ANM, its variables The distribution generation expression is as follows: (7) in, Representing variables The parent variable of For any function, different variables may correspond to different ; is independent Gaussian noise.

[0049] According to the different function forms and noise distributions, the model can be divided into linear Gaussian model, nonlinear Gaussian model, linear non-Gaussian model and nonlinear non-Gaussian model.

[0050] For the pairs of gene variables involved and In the simple case, based on the assumption that the cause and the noise are independent of each other, the causal relationship between the two can be tested by a nonlinear Gaussian model. In this embodiment, noise is estimated by Gaussian Process Regression (GPR), which is a classic non-parametric Bayesian regression method that satisfies the nonlinear relationship between the variable and the Gaussian noise, and can provide uncertainty measurement for the prediction results. The residual of the result obtained by GPR is regarded as noise. In this embodiment, the independence between the noise and the causal variable is verified by the Hilbert Schmidt Independence Criterion (HSIC). HSIC is a kernel-based independence criterion whose main purpose is to measure the distribution difference between two variables. Its principle is to define the Hilbert-Schmidt cross-covariance operator on the reproducing kernel Hilbert space, and then derive the statistic for measuring independence. The obtained test statistic can be used as a causal measure. The final ANM calculation is as shown in formula (8): (8) After completing the pairwise causal discovery based on ANM, two confidence matrices, namely the first confidence matrix, are obtained from the linear and nonlinear processes respectively. and the second confidence matrix , if the two processes are combined, more robust results may be obtained.

[0051] Here, this embodiment uses Borda counting voting method, which is a simple and effective merging method, is shown in formula (9). In this embodiment, the weight parameter Set to 0.5.

[0052] (9) The specific implementation process is shown in Algorithm 3:

[0053] Different models typically use different evaluation methods and metrics to assess the performance of GRN construction methods. This example considers the performance evaluation of GRN construction as a binary classification task. This example uses the commonly used metrics area under the receiver operating characteristic curve (AUROC) and area under the precision-recall curve (AUPR) to evaluate the performance of the GRN construction method.

[0054] In this example, the existence of a regulatory relationship derived from the model is considered a positive sample, while the absence of a regulatory relationship is considered a negative sample. The AUROC (Automatic Regression Operator) (AUROC) measures the model's ability to rank positive samples ahead of negative samples. Model performance is measured at different thresholds, plotting an ROC curve. The horizontal axis of the ROC curve represents the false positive rate (FPR), while the vertical axis represents the true positive rate (TPR). The relationship between the two reveals the robustness of the model's ability to identify positive samples. A value closer to 1 indicates better model performance. Furthermore, a key advantage of AUROC is that it does not directly depend on the distribution between classes, making it an effective metric even in imbalanced datasets. In contrast, the AUPR (Automatic Regression Operator) focuses on the precision and recall of positive samples, obtained by plotting a precision-recall curve. In imbalanced datasets, positive samples are typically far fewer than negative samples. In this case, even random guessing will result in a relatively small chance of misclassifying a positive sample, resulting in an artificially high recall rate. Therefore, the AUPR can be used as a more stringent metric that better reflects the model's ability to identify positive samples. The calculation formulas for TPR, FPR, precision and recall are as follows: (10) Among them, the real example Indicates the number of truly existing control edges that are correctly identified; false positives Indicates the number of non-real regulatory edges that are incorrectly predicted to exist; true counterexamples Indicates the number of correctly identified real non-existent control edges; false negative examples Indicates the number of incorrectly identified regulatory edges that do not actually exist.

[0055] The following experiments were conducted to verify the superiority of the solution of this embodiment: 1. Dataset

[0056] In this example, single-cell transcriptome data and GBM clinical datasets are used to construct gene regulatory network performance evaluation tasks and effect evaluation tasks for complex real biological datasets, respectively.

[0057] Single-cell transcriptome data is derived from the BEELINE benchmark library, including transcriptome data of human embryonic stem cells (hESC) and mouse mesenchymal cells (mDC). This embodiment collects three types of real networks: cell type-specific networks (Specific), nonspecific networks (NonSpecific), and STRING networks. Among them, cell type-specific networks and nonspecific networks are derived from the research of Pratapa et al., while STRING networks are constructed based on functional connectivity data in the STRING database. In order to generate a data set for gene regulatory network inference, this embodiment performs gene screening on each RNA sequencing data set, and extracts subsets containing 500 and 1000 genes with the largest variations, respectively. By crossing each real network with the gene subset, 12 data sets are finally generated. These data sets are named according to the single cell type, the real network type, and the number of genes screened out (500 or 1000). For example, the data set "mDC-Specific-500" represents a data set of the top 500 significantly expressed genes screened out based on the cell type-specific network of mouse mesenchymal cells.

[0058] In order to evaluate the effect of the framework of this embodiment on complex biological data sets, this embodiment uses a GBM clinical data set provided by a hospital. This data set integrates the RNA expression data of GBM-related genes in different patient samples from the TCGA, Rembrandt and CGGA-301 databases, as well as the corresponding patient clinical survival information. In the preprocessing stage, this embodiment screens out genes with available expression levels and removes genes whose expression values ​​are zero in more than 90% of the samples. The final number of genes retained were 16,176, 19,596 and 19,379, respectively, of which the number of transcription factors were 1,315, 1,424 and 1,420, respectively. A detailed description of the data sets involved in this embodiment is shown in Table 1.

[0059] Table 1: Dataset attributes

[0060] 2. Performance evaluation.

[0061] To evaluate the performance of the framework in this embodiment, we apply it to 12 datasets from the BEELINE benchmark library. Figure 2Comparisons are presented with prediction networks generated by 13 other widely used algorithms (including GENIE3, NARROMI, PIDC, GRNBOOST2, SPCO, PLSNET, CLR, BMA, catnet, ANM, Scribe, SINCERITIES, and GRINCD). The results of the comparison methods are based on the work of Jiang et al. This comparison focuses on four datasets related to human embryonic stem cells and mouse mesenchymal cells, using AUROC and AUPR as evaluation metrics. Results for the remaining datasets are shown in supplementary figures.

[0062] The results presented in Figures 4(a) and 4(b) show that the framework of this embodiment consistently outperforms the comparison algorithms in terms of both AUPR and AUROC metrics, demonstrating that the model of this embodiment possesses superior predictive capabilities. Overall, this model achieved the best performance across all four datasets. In terms of AUPR, this model outperformed the second-ranked method, NARROMI, by 28%–194.9%. In terms of AUROC, this model outperformed NARROMI by 19.5% on the hESC-Specific-500 dataset and by 9.31% on the hESC-CellType-1000 dataset. On the mDC-CellType-500 dataset, this model's AUROC was 42.43% higher than the worst-performing method, BMA; on the mDC-CellType-1000 dataset, it outperformed the worst-performing method, CATNET, by 32.37%. In addition, the results in the figure show that a single ANM shows poor results in both indicators on each data set. This is because ANM only relies on the information between paired nodes for structure discovery, ignoring the information of adjacent nodes, resulting in information loss and thus a decrease in the accuracy of the results.

[0063] 3. Ablation experiment.

[0064] To illustrate the necessity of adopting a dual-channel strategy combining linear and nonlinear methods within this framework, this example conducts ablation experiments. Specifically, by varying the channel combination threshold, the results are tailored to consider only linear or nonlinear channels. Table 2 shows the ablation experiment results for the four datasets considered.

[0065] Table 2: Ablation experiment results

[0066] The results show that under the same parameter settings and the same data set, the dual-channel merging of linear and nonlinear results can effectively improve model performance.

[0067] 4. Application and analysis in brain glioma scenarios.

[0068] This example uses a GBM clinical dataset provided by Shandong Qilu Hospital to validate the effectiveness of this framework for large-scale, high-dimensional, and complex single-cell transcriptome data (ScRNA-seq). Specifically, the validation was conducted from two perspectives: 1. Based on the constructed GRN, we analyzed the proportion of genes associated with TF knockout data, observing how the results compare to traditional correlation methods; 2. Based on clinical survival data, we performed a survival analysis based on expression patterns.

[0069] For 1, this embodiment uses the open source database KnockTF 1 Knockout data for 61 relevant TFs were downloaded from [1]. The classic correlation method, GENIE3, was selected as a control. The framework of this example and the selected control were applied to the three GBM datasets. Next, for the regulatory networks constructed using each method, a top-k=500 approach was used, selecting the top 500 TGs corresponding to each TF. The resulting mapping sets were compared with the TF knockout data. The resulting percentages of relevant TGs are shown in Figures 5(a)-5(c).

[0070] The results show that the framework of this embodiment achieves a higher percentage of results on all three datasets than GENIE3. This phenomenon indicates that when traditional correlation-based methods are used with large-scale, high-dimensional, and complex biological data, the TF-TF and TF-TG relationship pairs they obtain are not as comprehensive and accurate as those obtained by causal methods. This further demonstrates that the method of this embodiment is more effective than correlation methods in constructing true biological relationships.

[0071] Regarding 2, this example performed survival analysis on each TF based on clinical survival data. The surv_cutpoint function in the R package was used to classify TF expression data into optimally high and low expression, with the minimum ratio set to 0.3. As shown in Figures 6(a)-6(c), nine TFs were identified whose expression significantly impacted patient survival across all three datasets: ZIC3, TEAD2, PATZ1, PRKCQ, GRHL3, PYGO2, ZNF614, ZNF37A, and ZNF274.

[0072] In this example, we propose a framework that combines graph neural networks and causal discovery methods to infer causal gene regulatory networks. This framework leverages ClusterGCN's powerful global dependency capture capabilities to learn gene embedding representations and combines it with an improved additive noise model to analyze linear and nonlinear causal relationships, addressing the interpretability and accuracy shortcomings of traditional correlation-based methods.

[0073] Experimental validation using the BEELINE benchmark dataset and real-world glioblastoma clinical data showed that the proposed method significantly outperformed 13 mainstream methods across 12 datasets, improving the AUPR and AUROC metrics by 28%–194.9% and 9.31%–42.43%, respectively. Furthermore, the framework successfully identified nine transcription factors significantly associated with GBM patient survival, and wet experiments further validated the accuracy and biological significance of these regulatory relationships.

[0074] This paper proposes a framework that combines graph neural networks with causal discovery methods to construct a causal GRN. First, a gene co-expression network is constructed using the Pearson correlation coefficient and mutual information. Cluster-GCN is used to achieve efficient gene embedding representation learning to capture global dependencies. Second, an improved additive noise model (ANM) is used to analyze linear and nonlinear causal relationships, and a dual-channel integration strategy is used to enhance the robustness of the results. Experiments demonstrate that this framework significantly outperforms 13 mainstream methods on 12 datasets from the BEELINE benchmark, with improvements in AUPR and AUROC metrics of 28% to 194.9% and 9.31% to 42.43%, respectively. In real-world glioblastoma (GBM) clinical data, the model successfully identified nine transcription factors significantly associated with patient survival. Further wet experiments validated the accuracy of the regulatory relationships and the biological significance of the predicted TFs in GBM patient survival. By improving the accuracy and interpretability of causal discovery in high-dimensional heterogeneous transcriptome data, this paper provides a new paradigm for causal GRN inference in real biological data and a more accurate and reliable tool for analyzing complex biological systems.

[0075] Example 2 The purpose of this embodiment is to provide a gene regulatory network construction system, including: A construction module is configured to: determine the linear relationship and the nonlinear relationship between genes according to the gene expression matrix using the Pearson correlation coefficient method and mutual information, and respectively construct a first gene co-expression network and a second gene co-expression network; an embedding representation module configured to: based on the first gene co-expression network and the gene expression matrix, and the second gene co-expression network and the gene expression matrix, respectively use a graph convolutional network model to capture global dependencies and perform gene embedding representation learning to obtain a first gene embedding matrix and a second gene embedding matrix; An embedding representation module is configured to: calculate a causal relationship score based on the first gene embedding matrix and the second gene embedding matrix according to different types of regulatory relationships using an additive noise model to obtain a first confidence matrix and a second confidence matrix, and combine the first confidence matrix and the second confidence matrix to obtain a causal confidence matrix of the gene; The gene regulatory network module is configured to construct a gene regulatory network based on the causal confidence matrix of the gene and the gene expression matrix.

[0076] In further embodiments, there is also provided: An electronic device includes a memory and a processor, and computer instructions stored in the memory and executed by the processor. When the computer instructions are executed by the processor, the method described in Example 1 is performed. For the sake of brevity, no further details are given here.

[0077] It should be understood that in this embodiment, the processor may be a central processing unit (CPU), or may be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), off-the-shelf field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor, etc.

[0078] The memory may include a read-only memory and a random access memory, and provides instructions and data to the processor. A portion of the memory may also include a non-volatile random access memory. For example, the memory may also store information about the device type.

[0079] A computer-readable storage medium is used to store computer instructions, and when the computer instructions are executed by a processor, the method described in embodiment 1 is performed.

[0080] The method in Example 1 can be directly implemented as being executed by a hardware processor, or by a combination of hardware and software modules within the processor. The software module can be located in a storage medium well-established in the art, such as random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, or registers. The storage medium is located in the memory, and the processor reads the information in the memory and, in conjunction with its hardware, completes the steps of the above method. To avoid repetition, a detailed description is not given here.

[0081] Those skilled in the art will appreciate that the units and algorithm steps of the various examples described in conjunction with this embodiment can be implemented using electronic hardware or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professionals and technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present invention.

[0082] Although the above describes the specific embodiments of the present invention in conjunction with the accompanying drawings, it is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art on the basis of the technical solution of the present invention without any creative work are still within the scope of protection of the present invention.

Claims

1. A method for constructing a gene regulatory network, characterized in that: include: According to the gene expression matrix, the linear and nonlinear relationships between genes were determined using the Pearson correlation coefficient method and mutual information, and the first and second gene co-expression networks were constructed respectively; Based on the first gene co-expression network and the gene expression matrix, and the second gene co-expression network and the gene expression matrix, respectively, a graph convolutional network model is used to capture global dependencies for gene embedding representation learning to obtain a first gene embedding matrix and a second gene embedding matrix; Based on the first gene embedding matrix and the second gene embedding matrix, according to different types of regulatory relationships, a causal relationship score is calculated using an additive noise model to obtain a first confidence matrix and a second confidence matrix, and the first confidence matrix and the second confidence matrix are combined to obtain a causal confidence matrix of the gene; A gene regulatory network is constructed based on the causal confidence matrix of the gene and the gene expression matrix.

2. A method for constructing a gene regulatory network according to claim 1, characterized in that: According to the gene expression matrix, the linear and nonlinear relationships between genes were determined using the Pearson correlation coefficient method and mutual information, and the first and second gene co-expression networks were constructed respectively. Specifically: The linear and nonlinear relationships between genes in the gene expression matrix were calculated using the Pearson correlation coefficient method and mutual information, respectively, to construct the Pearson correlation coefficient matrix and the mutual information correlation coefficient matrix; According to the constructed Pearson correlation coefficient matrix and mutual information correlation coefficient matrix, and compared with the set threshold, the first gene co-expression network and the second gene co-expression network are determined respectively.

3. A method for constructing a gene regulatory network according to claim 1, characterized in that: Based on the first gene co-expression network and the gene expression matrix, or the second gene co-expression network and the gene expression matrix, a graph convolutional network model is used to capture global dependencies for gene embedding representation learning, thereby obtaining a first gene embedding matrix or a second gene embedding matrix, specifically: Normalize the gene expression matrix; Dividing a graph consisting of the first gene co-expression network and the normalized gene expression matrix, or a graph consisting of the second gene co-expression network and the normalized gene expression matrix, into a plurality of subgraphs; The graph convolutional network model is used to obtain the embedding representation of each subgraph, and all corresponding subgraph embeddings are merged to obtain the first gene embedding matrix or the second gene embedding matrix.

4. A method for constructing a gene regulatory network according to claim 3, characterized in that: The graph convolutional network model is used to obtain the embedding representation of each subgraph, and all corresponding subgraph embeddings are merged to obtain the first gene embedding matrix or the second gene embedding matrix, specifically: For each subgraph, extract the features corresponding to the subgraph node index set and the subgraph adjacency matrix; Calculating a degree matrix of the subgraph according to the subgraph adjacency matrix, and normalizing the subgraph adjacency matrix using the calculated degree matrix of the subgraph; Based on the normalized subgraph adjacency matrix and the features corresponding to the subgraph node index set, the subgraph embedding is calculated, and the subgraph embedding is mapped to the corresponding position of the original graph embedding matrix according to the node index to obtain the embedding representation of each subgraph. All corresponding subgraph embeddings are merged to obtain the first gene embedding matrix or the second gene embedding matrix.

5. A method for constructing a gene regulatory network according to claim 1, characterized in that: Based on the first gene embedding matrix or the second gene embedding matrix, according to different types of regulatory relationships, the causal relationship score is calculated using the additive noise model to obtain a first confidence matrix or a second confidence matrix, specifically: If the relationship pair belongs to the first gene embedding matrix or the second gene embedding matrix, for the relationship pair of transcription factors, two regressions in different directions are performed respectively and regression residuals are obtained respectively. The causal scores in the two directions are calculated by verifying the independence between the residuals and the independent variables; For the relationship between transcription factors and target genes, a regression in one direction is performed to obtain the regression residual, and the causal score in one direction is calculated by verifying the independence between the residual and the independent variable; The calculated causal scores are integrated to obtain a first confidence matrix or a second confidence matrix.

6. A method for constructing a gene regulatory network according to claim 5, characterized in that: The independence between noise and causal variables was verified by the Hilbert-Schmidt independence criterion.

7. A method for constructing a gene regulatory network according to claim 1, characterized in that: The gene regulatory network is a weighted directed bipartite graph ,in, refers to a collection of transcription factors, refers to the target gene set; yes There are directed edge sets in is the causal confidence matrix.

8. A gene regulatory network construction system, characterized in that: include: A construction module is configured to: determine the linear relationship and the nonlinear relationship between genes according to the gene expression matrix using the Pearson correlation coefficient method and mutual information, and respectively construct a first gene co-expression network and a second gene co-expression network; an embedding representation module configured to: based on the first gene co-expression network and the gene expression matrix, and the second gene co-expression network and the gene expression matrix, respectively use a graph convolutional network model to capture global dependencies and perform gene embedding representation learning to obtain a first gene embedding matrix and a second gene embedding matrix; An embedding representation module is configured to: calculate a causal relationship score based on the first gene embedding matrix and the second gene embedding matrix according to different types of regulatory relationships using an additive noise model to obtain a first confidence matrix and a second confidence matrix, and combine the first confidence matrix and the second confidence matrix to obtain a causal confidence matrix of the gene; The gene regulatory network module is configured to construct a gene regulatory network based on the causal confidence matrix of the gene and the gene expression matrix.

9. An electronic device, characterized in that: The method comprises a memory and a processor, and computer instructions stored in the memory and executed on the processor, wherein when the computer instructions are executed by the processor, the method according to any one of claims 1 to 7 is completed.

10. A computer-readable storage medium, characterized in that Used to store computer instructions, which, when executed by a processor, complete the method according to any one of claims 1 to 7.

Images