A causal gene regulatory network inference method based on directed graph attention network

Abstract

The present invention discloses a causal gene regulatory network inference method based on a directed graph attention network, which belongs to the field of graph representation learning and systems biology technology. The present invention uses prior regulatory information to construct an initial directed causal graph, and uses expression data as the initial features of the nodes in the graph. Then, two directed graph attention networks are used to learn the embedded representations of genes as source nodes and target nodes in the directed causal graph respectively, and a residual network is used to learn the characteristic information of the genes themselves implied in the expression data; finally, the probability of potential directed regulatory relationships is calculated based on the transcription factor embedding, target gene embedding and self-embedding of the genes. The present invention learns and models gene regulatory networks from the perspective of directed regulatory relationships, thereby improving the inference accuracy of potential directed regulatory relationships.

Description

Technical Field

[0001] The present invention relates to a causal gene regulatory network inference method based on a directed graph attention network, and belongs to the technical fields of graph representation learning, systems biology, etc. Background Art

[0002] Gene regulatory network is a network structure with transcription factors and target genes as nodes and directed causal regulatory relationships between them as directed edges. Many biological processes such as cell differentiation and cancer formation can be reflected in gene regulatory networks. Using gene expression data to accurately infer potential regulatory relationships is a long-term research topic in the field of systems biology.

[0003] The rapid development of single-cell RNA (scRNA-seq) sequencing technology has led to an exponential growth in single-cell gene expression data. Researchers have obtained a large amount of gene expression data and verified high-confidence prior regulatory relationship data. Therefore, it is necessary to develop computational methods that can use these data to discover potential causal regulatory relationships between genes.

[0004] Since gene regulatory networks can be viewed as directed causal graphs, many methods have emerged that use graph neural networks to learn gene representations and infer the probability of regulation between gene pairs. However, the currently proposed gene regulation inference methods based on graph neural networks are all based on undirected graphs to obtain gene representations, which results in the inability to fully capture the directed regulatory information provided by the prior regulatory network and to model the directed causal relationship of the gene regulatory network. This directed regulatory relationship is crucial because it determines the order of gene intervention and effect observation when conducting biological experimental verification. Summary of the invention

[0005] In view of the problem that the existing gene regulation inference methods based on graph neural networks cannot extract directed regulatory information in the prior regulatory network and cannot model the directed causal relationship of the gene regulatory network, the present invention provides a causal gene regulatory network inference method based on a directed graph attention network.

[0006] The present invention is implemented by the following technical scheme: a causal gene regulatory network inference method based on a directed graph attention network. This method regards the regulatory relationship between genes as a directed causal relationship, constructs an initial directed causal graph using prior regulatory information, and uses expression data as the initial features of the nodes in the graph. Then, two directed graph attention networks are used to learn the embedded representation of genes as source nodes (transcription factors) and target nodes (target genes) in the directed causal graph, respectively, and a residual network is used to learn the characteristic information of the genes themselves implied in the expression data; finally, the probability of potential directed regulatory relationships is calculated based on the transcription factor embedding, target gene embedding and self-embedding of the genes.

[0007] The specific steps are as follows (based on potential regulatory gene pairs For example):

[0008] Step 1: Use prior regulatory information to construct an initial directed regulatory network , and use the expression data as the initial features of the nodes in the graph: , , p is the initial gene feature dimension in the expression data;

[0009] Step 2: Use the encoder to learn different embedding representations of gene nodes;

[0010] Step 2.1, in the directed graph In the above, the identity of node a as the source node is defined as the transcription factor (TF) in the regulatory network, and its identity as the target node is defined as the target gene (TG) in the regulatory network. In this way, two kinds of neighborhoods of node a can be defined according to their different identities (TF, TG). One kind is , represents the first-order outgoing neighborhood of a, in which the node is a TF, and the other type is , represents the first-order incoming neighborhood of a, and the node is used as TG; then, DGAT is used to aggregate the neighborhood information of a through the attention mechanism, and each node has two identities in the graph , the attention score from b to a is:

[0011]

[0012] in For the connection operation, are the neighborhoods of node a with two different identities, s represents the transcription factor, t represents the target gene, is the weight matrix corresponding to the identity, is a single-layer feedforward neural network; finally, the TF embedding of a and TG embedded The output of the directed graph attention network DGAT is obtained as follows:

[0013]

[0014] in, , d is the output dimension of DGAT; this aggregates the outgoing neighborhood information of node a as TF and the incoming neighborhood information of TG; for the TF embedding and TG embedding of node b, the algorithm is the same as a;

[0015] Step 2.2: Use a residual network to obtain the characteristics of the gene itself hidden in the expression data. For gene a, its own embedding The calculation formula is as follows:

[0016]

[0017] in, , , n is the output dimension of the residual network, L is the number of residual blocks, is the output of the j-th layer residual block, is the j-th layer weight vector, represents the nonlinear transformation in the residual block; the calculation method of gene b’s own embedding is the same as that of a;

[0018] Step 3: Predict the regulation probability of potential directional regulatory gene pairs:

[0019] Step 3.1, determine whether gene a will have a causal regulatory effect on gene b through two items, including the consistency of TF embedding of a and self embedding of b, and the consistency of TG embedding of b and self embedding of a;

[0020] First, this method assumes that if gene a regulates gene b, then the TF of a is embedded in Embedded with self of b is close in the embedding space, through and The inner product between can measure the possibility of a regulating b. Similarly, this method assumes that if gene b is regulated by a, then the TG embedding of b Embedded with a's self Approximation in embedding space, through and The inner product of can measure the possibility that b is regulated by a. and , and They are in different embedding spaces, so the distance between them cannot be measured directly, so we need to use the alignment matrix , Come ( ) is mapped to ( ) space, a single-layer neural network is used for alignment, expressed as: ( ); then, through the inner product Get the possibility of gene a regulating gene b, and at the same time use the inner product Get the possibility that gene b is regulated by gene a;

[0021] After the above steps, we measured the To integrate the results of these two perspectives, we use the summation operation to obtain the probability of the directed causal regulation between a and b. , defined as follows:

[0022]

[0023] in, is the Sigmoid function, and are the coefficients of the linear combination, ;

[0024] Step 3.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:

[0025]

[0026] Among them, E is the real set of regulatory edge information, Indicates the existence of a known directed causal regulatory relationship , There is no known directed causal regulatory relationship ;

[0027] Step 3.3, Model performance evaluation:

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

[0029] In order to solve the problem of inferring directed causal regulatory relationships, the present invention introduces a directed graph attention network, and uses the source nodes and target nodes in the directed edges to distinguish the different identities of genes in the directed regulatory network, namely, transcription factors and target genes, and obtains different embedding representations of the corresponding identities. At the same time, the residual network is also used to obtain the embedding representation of its own characteristics; finally, the inference of potential directed regulatory relationships is completed based on the transcription factor embedding, target gene embedding and self-embedding of the gene.

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

[0031] The gene regulatory network inference method proposed in the present invention models the gene regulatory network as a directed causal graph, and learns different embedded representations of nodes in the graph through a directed graph attention network, effectively retaining the directed regulatory information in the prior network. At the same time, for potential directed regulatory relationship gene pairs, two items are used to measure their regulatory probability, namely: the consistency of the transcription factor embedding of the upstream regulatory gene and the self-embedding of the downstream target gene, and the consistency of the target gene embedding of the downstream target gene and the self-embedding of the upstream regulatory gene. It helps to clarify the direction of the regulatory relationship and improve the prediction performance of the potential regulatory relationship. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 A flow chart of a gene regulatory network inference method based on a directed graph attention network proposed in the present invention;

[0033] Figure 2 AUROC values ​​of each method on mouse and human scRNA-seq datasets;

[0034] Figure 3 AUPRC values ​​of each method on mouse and human scRNA-seq datasets. DETAILED DESCRIPTION

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

[0036] Example 1

[0037] like Figure 1 As shown, a gene regulatory network inference method based on a directed graph attention network, this embodiment first uses prior regulatory information to construct an initial directed causal regulatory network, and uses expression data as the initial features of the nodes. Then, two directed graph attention networks are used to learn the features of genes as source nodes (transcription factors) and target nodes (target genes) in the directed graph, and a residual network is used to learn the feature information of the genes themselves implied in the expression data. Finally, the encoder is used to calculate the consistency between the transcription factor embedding of the upstream gene and the self-embedding of the downstream gene in the potential directed regulatory gene pair, and the consistency between the target gene embedding of the downstream gene and the self-embedding of the upstream gene, and the two consistencies are combined to judge the potential directed regulation probability, so as to realize the inference of the directed gene regulatory network.

[0038] The specific steps include the following (using the potential regulatory gene pairs in the dataset For example):

[0039] Step 1: Use prior regulatory information to construct an initial directed regulatory network , and use the expression data as the initial features of the nodes in the graph: , , p is the initial gene feature dimension in the expression data.

[0040] Step 2: Use the encoder to learn different embedding representations of gene nodes.

[0041] Step 2.1, in the directed graph In the above, the identity of node a as the source node is defined as the transcription factor (TF) in the regulatory network, and its identity as the target node is defined as the target gene (TG) in the regulatory network. In this way, two kinds of neighborhoods can be defined according to the different identities of node a (TF, TG). One kind is , represents the first-order outgoing neighborhood of a, in which the node is a TF, and the other type is , represents the first-order incoming neighborhood of a, and the node is used as TG. Then, DGAT is used to aggregate the neighborhood information of a through the attention mechanism. Each node has two identities in the graph. , the attention score from b to a is:

[0042]

[0043] in For the connection operation, are the neighborhoods of node a with two different identities, s represents the transcription factor, t represents the target gene, is the weight matrix corresponding to the identity, is a single-layer feedforward neural network. Finally, the TF embedding of a and TG embedded It is obtained from the output of the directed graph attention network (DGAT) as follows:

[0044]

[0045] in, , d is the output dimension of DGAT. This aggregates the outgoing neighborhood information of node a as TF and the incoming neighborhood information of TG. For the TF embedding and TG embedding of node b, the algorithm is the same as a.

[0046] Step 2.2: Use a residual network to obtain the characteristics of the gene itself hidden in the expression data. For gene a, its own embedding The calculation formula is as follows:

[0047]

[0048] in, , , n is the output dimension of the residual network, L is the number of residual blocks, is the output of the j-th layer residual block, The jth layer is the weight vector, represents the nonlinear transformation in the residual block. The embedding of gene b itself is calculated in the same way as a.

[0049] Step 3: Predict the regulatory probability of potential directed regulatory gene pairs.

[0050] Step 3.1. Use two items to determine whether gene a will have a causal regulatory effect on gene b, including the consistency of a's TF embedding and b's self embedding, and the consistency of b's TG embedding and a's self embedding.

[0051] First, this method assumes that if gene a regulates gene b, then the TF of a is embedded in Embedded with self of b is close in the embedding space, through and The inner product between can measure the possibility of a regulating b. Similarly, this method assumes that if gene b is regulated by a, then the TG embedding of b Embedded with a's self Approximation in embedding space, through and The inner product of can measure the possibility that b is regulated by a. and , and They are in different embedding spaces, so the distance between them cannot be measured directly, so we need to use the matrix , Come ( ) is mapped to ( ), here we use a single-layer neural network for alignment, which can be expressed as: ( ). After that, we can use the inner product Get the possibility of gene a regulating gene b, and at the same time use the inner product Get the possibility that gene b is regulated by a.

[0052] After the above steps, we measured the To integrate the results of these two perspectives, we use the summation operation to obtain the probability of the directed causal regulation between a and b. , defined as follows:

[0053]

[0054] in, is the Sigmoid function, and are the coefficients of the linear combination, .

[0055] Step 3.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:

[0056]

[0057] Among them, E is the real set of regulatory edge information, Indicates the existence of a known directed causal regulatory relationship , There is no known directed causal regulatory relationship ;

[0058] Step 3.3, model performance evaluation;

[0059] Step 3.3.1. Evaluation indicators: The area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC) were used as evaluation indicators.

[0060] Step 3.3.2, Experimental data set: scRNA-seq data sets of seven cell lines from humans and mice were selected to evaluate the performance of the model, 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 data set, all transcription factors with corrected p-values ​​less than 0.01 and the top (500 / 1000) significantly changed target genes were selected for regulatory relationship inference. At the same time, three real networks were selected as prior regulatory information, which came from the functional interaction networks recorded in the STRING database, non-specific ChIP-seq, and cell type-specific ChIP-seq. The detailed information is shown in Table 1:

[0061] Table 1 Dataset information statistics

[0062]

[0063] Step 3.3.3 Experimental results: For the above data sets, the positive and negative samples of the training set, validation set and test set are divided according to the network density; AUROC and AUPRC are selected as evaluation indicators, and 5 experiments are conducted on each data set and the results are finally averaged. The results are as follows: Figure 2 and Figure 3 As shown, DGATGRN is the method proposed by the present invention.

[0064] The specific implementation modes of the present invention are described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above implementation modes, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A causal gene regulatory network inference method based on directed graph attention network, characterized by: The initial directed causal graph is constructed using prior regulatory information, and the expression data is used as the initial features of the nodes in the graph. Then, two directed graph attention networks are used to learn the embedded representations of genes as source nodes and target nodes in the directed causal graph, respectively. At the same time, the residual network is used to learn the characteristic information of the genes themselves implied in the expression data. Finally, the probability of potential directed regulatory relationships is calculated based on the transcription factor embedding, target gene embedding and self-embedding of the gene. The specific steps are as follows: Step 1: Use prior regulatory information to construct an initial directed regulatory network , and use the expression data as the initial features of the nodes in the graph: , , p is the initial gene feature dimension in the expression data; Step 2: Use the encoder to learn different embedding representations of gene nodes; Step 2.1, in the directed graph In the above example, the identity of node a as a source node is defined as a transcription factor in the regulatory network, and its identity as a target node is defined as a target gene in the regulatory network. According to the different identities of node a, two types of neighborhoods are defined: , represents the first-order outgoing neighborhood of a, in which the node is a TF, and the other type is , represents the first-order incoming neighborhood of a, and the node is used as TG; then, DGAT is used to aggregate the neighborhood information of a through the attention mechanism, and each node has two identities in the graph , the attention score from b to a is: ； in For the connection operation, are the neighborhoods of node a with two different identities, s represents the transcription factor, t represents the target gene, is the weight matrix corresponding to the identity, is a single-layer feedforward neural network; finally, the TF embedding of a and TG embedded The output of the directed graph attention network DGAT is obtained as follows: ； in, , d is the output dimension of DGAT; this aggregates the outgoing neighborhood information of node a as TF and the incoming neighborhood information of TG; for the TF embedding and TG embedding of node b, the algorithm is the same as a; Step 2.2: Use a residual network to obtain the characteristics of the gene itself hidden in the expression data. For gene a, its own embedding The calculation formula is as follows: ； in, , , n is the output dimension of the residual network, L is the number of residual blocks, is the output of the j-th layer residual block, is the j-th layer weight vector, represents the nonlinear transformation in the residual block; the calculation method of gene b’s own embedding is the same as that of a; Step 3: Predict the regulation probability of potential directional regulatory gene pairs: Step 3.1, determine whether gene a will have a causal regulatory effect on gene b through two items, including the consistency of TF embedding of a and self embedding of b, and the consistency of TG embedding of b and self embedding of a; First, by aligning the matrix , Come ( ) is mapped to ( ) space, a single-layer neural network is used for alignment, expressed as: ( ); then, through the inner product Get the possibility of gene a regulating gene b, and at the same time use the inner product Get the possibility that gene b is regulated by a; After the above steps, we measured the The probability of a directed causal regulation between a and b is obtained by summing the two. , defined as follows: ； in, is the Sigmoid function, and are the coefficients of the linear combination, ; Step 3.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: ； Among them, E is the real set of regulatory edge information, Indicates the existence of a known directed causal regulatory relationship , There is no known directed causal regulatory relationship ; Step 3.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