Gene regulatory network inference method based on second-order multi-channel directed graph convolution

Abstract

The invention provides a gene regulation network inference method based on second-order multi-channel directed graph convolution, and belongs to the technical fields of graph representation learning, system biology and the like. Firstly, a first-order proximity matrix, a second-order in-degree proximity matrix and a second-order out-degree proximity matrix are constructed respectively, an existing method is expanded through the second-order proximity matrix, secondly, feature information implied in the three matrixes and expression data is extracted by using multiple graph convolution networks respectively, and the feature information of the three matrixes and the expression data is extracted. And then fusion operation is carried out to fuse features extracted by the multi-channel directed graph convolutional network. And finally, respectively extracting embedding expressions of transcription factors and target genes in potential regulatory gene pairs by using two channels, and obtaining a predicted regulatory probability through dot product operation. According to the method, the prior regulation and control information and the directional regulation and control information implied in the expression data are extracted, and the perception domain of the graph convolution is expanded from the first order to the second order, so that richer feature information can be extracted, and the inference accuracy of the potential directional regulation and control relationship is improved.

Description

Technical Field

[0001] The present invention relates to a method for inferring gene regulatory networks based on second-order multi-channel directed graph convolution, and belongs to the technical fields of graph representation learning, systems biology, etc. Background Art

[0002] Gene regulatory networks are important tools for describing the regulatory relationships between genes, and are widely used in revealing gene expression regulatory mechanisms, exploring cell differentiation, and the formation, prevention, and treatment of diseases. By analyzing the complex interaction patterns between genes through gene regulatory networks and mining key genes related to specific functions, it is possible to save time and funds for biological researchers and improve efficiency.

[0003] The rapid development of single-cell RNA (scRNA-seq) sequencing technology has led to an exponential growth in single-cell gene expression data, and researchers have obtained a large number of gene expression datasets. Therefore, there is a need to develop computational methods that can use these data to discover potential regulatory relationships between genes. Currently, various regulatory network inference methods based on gene expression data have been proposed. Since gene regulatory networks can be naturally constructed as a graph structure, many methods have emerged that use graph neural networks to learn gene representations and infer the regulatory probabilities between gene pairs. However, the currently proposed gene regulatory inference methods based on graph neural networks all obtain gene representations based on undirected graphs, and cannot extract directed regulatory information. At the same time, it has also been found that due to limited known regulatory information, there may not be known direct connections between many genes with the same or similar expression patterns, and current methods usually only extract information from the first-order neighborhood, which may lead to the inability to fully capture more abundant information in the prior network and expression data. Summary of the Invention

[0004] Aiming at the problems that the existing gene regulatory inference methods based on graph neural networks cannot extract the directed regulatory information in the prior regulatory network, and at the same time, there is also the problem of only extracting the first-order neighborhood and limited feature extraction ability in the case of insufficient prior regulatory information, the present invention provides a method for inferring gene regulatory networks based on second-order multi-channel directed graph convolution.

[0005] The present invention is realized through the following technical solutions: A method for inferring gene regulatory networks based on second-order multi-channel directed graph convolution. Based on the prior regulatory network, this method first constructs three matrices, namely the first-order proximity matrix, the second-order in-degree proximity matrix, and the second-order out-degree proximity matrix, and uses the second-order proximity matrix to expand the existing method, which not only retains the directional characteristics of the directed graph, but also can expand the receptive field of graph convolution. Then, the local and global structures of the directed graph are more accurately characterized through the directional Laplacian matrix, and the feature representations of genes are obtained and used for regulatory relationship inference, further improving the performance of network inference and the modeling ability for complex regulatory patterns.

[0006] The specific steps are as follows:

[0007] Step1. Use the initial directed regulatory network provided by the prior regulatory information to construct the first-order and second-order proximity matrices;

[0008] Step1.1 First-order proximity: First-order proximity refers to the local pairwise structure between vertices in a graph. To model the first-order proximity between vertices, define the first-order proximity degree between vertex and vertex in each edge (i, j) of graph G as follows:

[0009] (1)

[0010] where is the symmetric matrix representing the adjacency matrix of the original graph G. If there is no edge from vertex to vertex or from vertex to vertex , then , is the initial adjacency matrix, and d is the number of genes. The restriction on the directed Figure 1 first-order proximity here is relatively loose, and the symmetric matrix is used to replace the original matrix, which inevitably loses some directed information. For this part of the missing information, another method will be used to retain it, that is, second-order proximity.

[0011] Step1.2 Second-order proximity. Second-order proximity assumes that if two vertices have multiple common neighbors, they are similar. In an actual regulatory network, since the cost of obtaining known regulatory relationships is relatively high, there are no direct edges between many very similar genes. Therefore, it is necessary to establish a second-order proximity matrix to connect similar nodes. The second-order proximity degree between vertices is determined by the sum of the normalized weights of the edges connected to the nodes in their shared neighborhood. In a directed graph G, for vertex and vertex , their second-order in-degree proximity and second-order out-degree proximity are defined as follows:

[0012] (2)

[0013] (3)

[0014] where is the adjacency matrix representing the original graph G, , V is the set of vertices in the original graph G.

[0015] The larger it is, the higher the second-order in-degree similarity between two vertices. Similarly, The larger it is, the higher the second-order out-degree similarity between two vertices. If there are no common vertices between vertex and vertex , then their second-order proximity is 0. In a directed graph, the edges adjacent to vertex and vertex and vertex are in pairs. Therefore, , . Although and are symmetric, the directed information in the original graph G is retained through their construction process.

[0016] Step 2: Use multi-channel directed graph convolution to obtain different embedding representations of gene nodes.

[0017] The formula for ordinary graph convolution is as follows:

[0018] (4)

[0019] where is the feature matrix, d is the number of genes, p is the initial feature dimension of the genes, is the trainable parameter matrix, h is the output dimension, , is the adjacency matrix, is the identity matrix, is 's degree matrix, where , and H is the output result. This kind of graph convolution is only applicable to undirected graphs. When dealing with directed graphs, it can only transform the directed graph into an undirected graph first and construct a symmetric Laplacian matrix, which will undoubtedly lead to the loss of some directed information.

[0020] In the previous step, the first-order proximity and second-order proximity were defined, and three proximity matrices were also obtained, namely the first-order proximity matrix, the second-order in-degree proximity matrix, and the second-order out-degree proximity matrix. In formula (4), the adjacency matrix A of the graph stores the information of the graph and provides a receptive field for to realize the transformation from the original feature X to Z. Since the three matrices obtained in the previous steps are similar to A and are all symmetric matrices, the first-order proximity matrix provides a receptive field similar to 1-hop, and the second-order in-degree proximity matrix and the second-order out-degree proximity matrix provide receptive fields similar to 2-hop, while retaining the directed information in the original graph. Therefore, based on ordinary graph convolution, this step defines the first-order proximity convolution Second-order in-degree proximity convolution and second-order out-degree proximity convolution are as follows:

[0021] (5)

[0022] (6)

[0023] (7)

[0024] Among them, , and are the outputs of three different convolution operations respectively. Through the above multi-channel proximity convolution operations, , and not only capture the features of rich first-order and second-order neighbors, and but also retain the directed structure information in the original graph, thus realizing the multi-channel directed graph convolution operation.

[0025] Step 3: Fusion operation. The above steps realize the directed graph convolution operation through multi-channel graph convolution operations and also obtain three different features. Next, they need to be combined to obtain the fused gene feature representation , and the specific formula is as follows:

[0026] (8)

[0027] Among them, represents the fusion function, which can be customized according to requirements and can be a summation function, a normalization function, a concatenation function, or a fusion function based on attention weights. The fusion operation function specifically used in the present invention is as follows:

[0028] (9)

[0029] Among them, Concat represents matrix concatenation operation, and are the weight parameters that control the output features of the second-order in-degree proximity convolution and the output features of the second-order out-degree proximity convolution, and can be freely combined according to the characteristics of the network topology. For example, in a graph with fewer second-order neighbors, the values of and can be appropriately reduced, and more first-order information can be used. Similarly, in the second-order neighbors, if there are fewer in-degree neighbors and more out-degree neighbors, then by adjusting and Adjust the amplitude of their impact on the whole according to the ratio, and vice versa. After the above steps, the comprehensive gene feature representation Z that integrates first-order and second-order information as well as directed structural information can be obtained.

[0030] Step4. Prediction of potential directed regulatory relationships.

[0031] Step4.1. After obtaining the comprehensive gene feature representation Z, it is necessary to predict the regulatory relationship between genes based on this. For a pair of potential regulatory gene pairs (i, j), the vector representations corresponding to them in Z are input into two identical channels in this step, which are used to obtain the low-dimensional transcription factor embedding and target gene embedding respectively. Each channel is composed of fully connected layers, and the specific representation is as follows:

[0032] (10)

[0033] (11)

[0034] Among them, LeakyReLU is the activation function, 、 、 and are the transposes of four weight matrices respectively, 、 、 and are four bias vectors respectively, is the fused feature representation of gene i, is the fused feature representation of gene j, is the final embedding representation vector of gene i, is the final embedding representation vector of gene j, and n is the dimension of the final embedding representation of the gene; finally, the directed regulatory probability between gene i and gene j is predicted through dot product , and the specific is as follows:

[0035] (12)

[0036] Step4.2. Take the existing interaction pairs as positive samples, and randomly select the non-existing interaction pairs as negative samples. Use the Adam optimizer to train the model, and the binary cross-entropy is used as the loss function. The specific is as follows:

[0037] (13)

[0038] Among them, M is the number of transcription factor-target gene pairs, i represents the i-th pair of transcription factor-target gene pairs, represents the label of the i-th pair of transcription factor-target gene pairs, represents the prediction result of the i-th pair of transcription factor-target gene pairs.

[0039] Step 4.3. Model performance evaluation;

[0040] Evaluation metrics: The area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC) are used as evaluation metrics.

[0041] To solve the problems of extracting directed regulatory information and the limited receptive field and insufficient information capture ability caused by only aggregating first-order neighborhoods due to limited prior regulatory information, the present invention introduces second-order multi-channel directed graph convolution. It not only extracts first-order neighborhood information but also expands the receptive field of the graph convolutional network by constructing a second-order proximity matrix, expands the scope of information aggregation to the second order, then obtains comprehensive gene representations through a fusion operation, and finally obtains the embedding representations of potential directed regulatory gene pairs for upstream transcription factors and downstream target genes through two channels respectively, and uses the dot product to predict the regulatory probability.

[0042] The beneficial effects of the present invention are:

[0043] The gene regulatory network inference method proposed by the present invention not only retains the directed regulatory information in the prior network by constructing first-order and second-order in-degree and out-degree proximity matrices, but also improves the receptive field of the graph convolutional network from the first order to the second order, enabling the convolutional network to extract more feature information hidden in the prior network and expression data. In the case of limited prior regulatory information, this method helps to extract richer information and clarify the direction of regulatory relationships, improving the prediction performance of potential regulatory relationships. Brief Description of the Drawings

[0044] Figure 1 It is a flowchart of a gene regulatory network inference method based on second-order multi-channel directed graph convolution proposed by the present invention;

[0045] Figure 2 It is the AUROC values of each method on mouse and human scRNA-seq datasets;

[0046] Figure 3 It is the AUPRC values of each method on mouse and human scRNA-seq datasets;

[0047] Figure 2 and Figure 3 In, sMDGCGRN is the method adopted by the present invention;

[0048] Others are methods adopted in the prior art:

[0049] scMGATGRN is a gene regulatory inference method for multi-view graph attention networks; GENELink is a method for predicting potential regulatory relationships based on graph attention networks; GNE is a method for encoding gene expression profiles and network topologies based on MLP to predict gene dependencies; DeepSEM is the first method to apply deep learning to gene regulatory network prediction; GRNBoost2 is an efficient gene regulatory network inference algorithm based on gradient boosting machines; GENIE3 is an inference method based on random forests. Detailed implementation manners

[0050] The present invention will be further described below in conjunction with embodiments.

[0051] Embodiment 1

[0052] As Figure 1 shown, a gene regulatory network inference method based on second-order multi-channel directed graph convolution. In this embodiment, first, the first-order proximity matrix, the second-order in-degree proximity matrix, and the second-order out-degree proximity matrix are respectively constructed using prior regulatory information. The existing methods are extended through the second-order proximity matrix, which can not only retain the directional features of the directed graph but also expand the receptive field of graph convolution, thereby extracting more features. Secondly, multiple graph convolution networks are used to extract the feature information hidden in these three matrices respectively. Then, a fusion operation is performed to fuse the features extracted by the multi-channel directed graph convolution network. Finally, two channels are used to extract the embedding representations of the transcription factor and the target gene in the potential regulatory gene pairs respectively, and the predicted regulatory probability is obtained through dot product operation.

[0053] Specifically, it includes the following steps:

[0054] Step1. Use the initial directed regulatory network provided by prior regulatory information , and construct the first-order and second-order proximity matrices.

[0055] Step1.1. First-order proximity. First-order proximity refers to the local pairwise structure between vertices in a graph. To model the first-order proximity between vertices, define the first-order proximity degree between vertex and vertex in each edge (i, j) in graph G as follows:

[0056] (1)

[0057] where is the symmetric matrix representing the adjacency matrix of the original graph G. If there is no edge from vertex to vertex or from vertex to vertex , then , is the initial adjacency matrix, and d is the number of genes. Here, the restriction on the Figure 1 second-order proximity is relatively loose. A symmetric matrix is used to replace the original matrix, which inevitably loses some directed information. For this part of the missing information, another method will be used to retain it, that is, the second-order proximity.

[0058] Step1.2 Second-order proximity. The second-order proximity assumes that if two vertices have multiple common neighbors, they are similar. In an actual regulatory network, since the cost of obtaining known regulatory relationships is relatively high, there are no direct edges between many very similar genes. Therefore, it is necessary to establish a second-order proximity matrix to connect similar nodes. The second-order proximity degree between vertices is determined by the sum of the normalized weights of the edges connected to their shared neighborhood nodes. In a directed graph G, for vertices and vertex , their second-order in-degree proximity and second-order out-degree proximity are defined as follows:

[0059] (2)

[0060] (3)

[0061] where represents the adjacency matrix of the original graph G, , and V is the set of vertices in the original graph G.

[0062] The larger is, the higher the second-order in-degree similarity between the two vertices. Similarly, and vertex The larger is, the higher the second-order out-degree similarity between the two vertices. If there are no common vertices between vertex and vertex , their second-order proximity degree is 0. In a directed graph, the edges adjacent to vertex , , although and are symmetric, the directed information in the original graph G is retained through their construction process.

[0063] Step2. Obtain different embedding representations of gene nodes using multi-channel directed graph convolution.

[0064] The general graph convolution formula is as follows:

[0065] (4)

[0066] Among them, is the feature matrix, d is the number of genes, and p is the initial feature dimension of the genes. is the trainable parameter matrix, h is the output dimension. , is the adjacency matrix. is the identity matrix. is 's degree matrix, where , and H is the output result. This kind of graph convolution is only applicable to undirected graphs. When dealing with directed graphs, it can only first convert the directed graph into an undirected graph and construct a symmetric Laplacian matrix, which will undoubtedly result in the loss of some directed information.

[0067] In the previous step, first-order proximity and second-order proximity were defined, and three proximity matrices were also obtained, namely the first-order proximity matrix, the second-order in-degree proximity matrix, and the second-order out-degree proximity matrix. In formula (4), the adjacency matrix A of the graph stores the information of the graph and provides a receptive field for , thus realizing the transformation from the original feature X to Z. Since the three matrices obtained in the previous steps are similar to A and are all symmetric matrices, the first-order proximity matrix provides a receptive field similar to 1-hop, and the second-order in-degree proximity matrix and the second-order out-degree proximity matrix provide receptive fields similar to 2-hop, while retaining the directed information in the original graph. Therefore, based on ordinary graph convolution, this chapter defines the first-order proximity convolution , the second-order in-degree proximity convolution , and the second-order out-degree proximity convolution as follows:

[0068] (5)

[0069] (6)

[0070] (7)

[0071] Among them, , , and are the outputs of three different convolution operations respectively. Through the above multi-channel proximity convolution operations, , , and not only capture the features of rich first-order and second-order neighbors, , and also retain the directed structure information in the original graph, thereby realizing the multi-channel directed graph convolution operation.

[0072] Step 3. Fusion operation. The above steps implement the directed graph convolution operation through the multi-channel graph convolution operation, and three different features are also obtained. Next, they need to be combined to obtain the fused gene feature representation. , and the specific formula is as follows:

[0073] (8)

[0074] Among them, represents the fusion function, which can be customized according to needs. It can be a summation function, a normalization function, a concatenation function, or a fusion function based on attention weights. The specific fusion operation function used in this method is as follows:

[0075] (9)

[0076] Among them, Concat represents the matrix concatenation operation, and are the weight parameters that control the output features of the second-order in-degree neighbor convolution and the second-order out-degree neighbor convolution, and can be freely combined according to the characteristics of the network topology. For example, in a graph with fewer second-order neighbors, the values of and can be appropriately reduced, and more first-order information can be used. Similarly, in the second-order neighbors, if there are fewer in-degree neighbors and more out-degree neighbors, the ratio of and can be adjusted to adjust their influence on the whole, and vice versa. After the above steps, the comprehensive gene feature representation Z that integrates the first-order, second-order information, and the directed structure information can be obtained.

[0077] Step 4. Prediction of potential directed regulatory relationships.

[0078] Step 4.1. After obtaining the comprehensive gene feature representation Z, it is necessary to predict the regulatory relationship between genes based on this. For a pair of potential regulatory gene pairs (i, j), in this chapter, their corresponding vector representations in Z are input into two identical channels to obtain the low-dimensional transcription factor embedding and target gene embedding respectively. Each channel is composed of a fully connected layer, and the specific representation is as follows:

[0079] (10)

[0080] (11)

[0081] Among them, LeakyReLU is the activation function, , , and are the transposes of four weight matrices respectively. , , and are four bias vectors respectively. is the fusion feature representation of gene i. is the fusion feature representation of gene j. is the final embedded representation vector of gene i. is the final embedded representation vector of gene j, and n is the dimension of the final embedded representation of the gene. Finally, the directed regulation probability between gene i and gene j is predicted through dot product as follows:

[0082] (12)

[0083] Step4.2. Take the existing interaction relation pairs as positive samples, and randomly select the non-existing interaction relation pairs as negative samples. Use the Adam optimizer to train the model, and binary cross-entropy as the loss function, as follows:

[0084] (13)

[0085] where M is the number of transcription factor-target gene pairs, and i represents the i-th pair of transcription factor-target gene pairs. represents the label of the i-th pair of transcription factor-target gene pairs. represents the prediction result of the i-th pair of transcription factor-target gene pairs.

[0086] Step4.3. Model performance evaluation;

[0087] Step4.3.1. Evaluation metrics: The area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC) are used as evaluation metrics.

[0088] Step4.3.2. Experimental Datasets: The scRNA-seq datasets of seven cell lines of humans and mice were selected to evaluate the model performance, namely: mouse embryonic stem cells (mESC), mouse dendritic cells (mDC), mouse erythroid hematopoietic stem cells (mHSC-E), mouse hematopoietic stem cells with granulocyte-monocyte lineage (mHSC-GM), mouse hematopoietic stem cells with lymphoid lineage (mHSC-L), human embryonic stem cells (hESC), and human mature hepatocytes (hHEP). For each dataset, all transcription factors with corrected p-values less than 0.01 and the top (500 / 1000) significantly changed target genes (TFs+500 / TFs+1000) were selected for regulatory relationship inference. At the same time, the functional interaction network recorded in the non-cell-specific ChIP-seq database (Non-specific CHIP-seq) was selected as prior information. The detailed information is shown in Table 1:

[0089] Table 1 Dataset statistical information corresponding to the non-cell-specific ChIP-seq database

[0090] Step4.3.3. Experimental Results: For the above datasets, the positive and negative sample ratios were balanced in the training set and the validation set. At the same time, the positive and negative samples of the test set were divided according to the network density; AUROC and AUPRC were selected as evaluation indicators, and 5 experiments were conducted on each dataset, and the results were finally averaged. The results are as Figure 2 and Figure 3 shown, where sMDGCGRN is the method proposed in the present invention.

[0091] The specific embodiments of the present invention have been described in detail above in conjunction with the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those of ordinary skill in the art, various changes can be made without departing from the spirit of the present invention.

Claims

1. A gene regulatory network inference method based on second-order multi-channel directed graph convolution, characterized by: A first-order proximity matrix representing the direct regulatory relationship is constructed using prior regulatory information. At the same time, a second-order out-degree proximity matrix and a second-order in-degree proximity matrix are constructed based on the known directed regulatory similarities between genes. These three matrices together with the expression data are then input into a multi-channel graph convolutional network to extract different feature information. The obtained multi-channel feature information is fused to obtain a comprehensive gene representation. Subsequently, two channels consisting of fully connected layers are used to obtain the embedded representations of transcription factors and target genes in potential regulatory gene pairs. Finally, a dot product operation is performed to obtain the potential directed regulation probability.

2. A gene regulatory network inference method based on second-order multi-channel directed graph convolution according to claim 1, characterized in that The specific steps are as follows: Step 1: Use the initial directed regulatory network provided by prior regulatory information , construct first-order and second-order proximity matrices; Step 1.1, first-order proximity: define the vertex of each edge (i, j) in graph G and vertices The first-order proximity between As shown below: （1）； in is the adjacency matrix representing the original graph G If there is no symmetric matrix from the vertex To the top Or from the vertex To the top The edge of ; is the initial adjacency matrix, d is the number of genes; Step 1.2, Second-order proximity: In the directed graph G, for the vertex and vertices , whose second-order in-degree is adjacent to and second-order outdegree proximity The definition is as follows: （2）； （3）； in, is the adjacency matrix representing the original graph G, , V is the vertex set in the original graph G; Step 2: Use multi-channel directed graph convolution to obtain different embedding representations of gene nodes: define first-order neighbor convolution , second-order in-degree neighbor convolution and second-order out-degree neighbor convolution , as shown below: （5）； （6）； （7）； in, , and They are the outputs of three different channels; is the feature matrix, d is the number of genes, p is the initial feature dimension of the gene, is a trainable parameter matrix, h is the output dimension, , is the adjacency matrix, is the identity matrix, ; Step 3: Fusion operation: （9）； Among them, Concat represents the matrix concatenation operation, and It is the weight parameter that controls the output features of the second-order in-degree neighboring convolution and the second-order out-degree neighboring convolution. After the above steps, the comprehensive gene feature representation Z that integrates the first-order, second-order information, and directed structural information is obtained; Step 4: Prediction of potential directional regulatory relationships: Step 4.

1. For a pair of potential regulatory gene pairs (i, j), input their corresponding vector representations in Z into two identical channels, as follows: （10）； （11）； Among them, LeakyReLU is the activation function, , , and are the transposes of the four weight matrices, , , and They are four bias vectors, is the fusion feature representation of gene i, is the fusion feature representation of gene j, is the final embedding representation vector of gene i, is the final embedding representation vector of gene j, and n is the final embedding representation dimension of the gene; finally, the dot product is used to predict the directional regulation probability between gene i and gene j , as follows: （12）； Step 4.2: Take the existing interaction relationship pairs as positive samples, and randomly select the non-existent interaction relationship pairs as negative samples. Use the Adam optimizer to train the model, and use the binary cross entropy as the loss function. The details are as follows: （13）； Where M is the number of transcription factor-target gene pairs, i represents the i-th transcription factor-target gene pair, represents the label of the i-th transcription factor-target gene pair, represents the prediction result of the i-th transcription factor-target gene pair; Step 4.3, model performance evaluation; Evaluation indicators: The area under the receiver operating characteristic curve and the area under the precision-recall curve were used as evaluation indicators.

Images