A gene regulatory network topology identification method and system based on random variation Bayesian and a storage medium

Abstract

The present application relates to the technical field of biological gene expression, and particularly relates to a gene regulation network topology identification method and system based on random variation Bayesian and a storage medium. The steps are as follows: 1) obtaining a gene expression dataset to be identified; 2) reconstructing a gene regulation network model by using a dynamic structure function DSF based on a state space model; 3) estimating parameters in the gene regulation network model by using a method based on random variation Bayesian; 4) identifying the topology of the gene regulation network by using a forward selection method, updating the network model by using an ARD variable, selecting a model structure by using a lower bound function J, and drawing a gene regulation network topology graph. The present application updates the posterior distribution of a global variable by using a natural gradient based on partial data, and the calculation cost is much lower than that of a classical VI method, so the present application is suitable for identifying the topology of a gene regulation network.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of biological gene expression technology, and particularly to a gene regulatory network topology identification method and system based on random variation Bayesian and a storage medium. BACKGROUND

[0002] Bioinformatics research shows that some specific physiological activities of some organisms are related to the expression profile values of some genes in the cells of the organisms, the expression value of a gene of an organism is affected by the expression values of other genes, and the expression value of the gene also affects the expression values of some other genes, and the regulatory relationship between the genes that restrict each other constitutes a complex gene regulatory network (GRN) in the field of bioinformatics. Identifying the gene regulatory network is the basis for further understanding the complex regulatory relationship between genes.

[0003] Using mathematical modeling methods and systems biology knowledge, analyzing gene expression data to reconstruct gene regulatory networks is a commonly used method for constructing gene regulatory networks at present, as shown in Figure 4 . Figure 4 In the figure, (a) is gene expression data, the row vector of the matrix represents genes (Gene, G), and the column vector represents samples (Sample, S); (b) is a graph of the regulatory relationship between genes constructed by using a network inference method.

[0004] At present, there are many gene regulatory network inference methods and models that have been developed, which can be roughly divided into correlation-based, Boolean network, differential equation model, and Bayesian network-based gene regulatory network inference methods.

[0005] Among them, Bayesian networks are widely used. Using Bayesian methods to predict the topology of gene regulatory networks has the following advantages:

[0006] 1. Bayesian networks are directed acyclic graphs that can reflect the dependency relationship between genes.

[0007] 2. Bayesian networks have certain processing capabilities for data noise.

[0008] 3. Bayesian networks are clear in expression, for example, if there is no connection between nodes, it means that there is no probability relationship between the nodes, indicating that there is no regulatory effect between the genes.

[0009] The purpose of using Bayesian network methods to identify gene regulatory networks is to find the probability relationship of the mutual regulation of genes and to determine whether the current gene in the network is regulated by other genes or inputs. The accuracy of the final identification result is reflected by the true positive rate.

[0010] Compared with other models, Bayesian network can find the optimal structure of the predicted network matching the real gene regulatory network, but traditional Bayesian network also has certain shortcomings, and due to its high computational time complexity, can only be applied to small-scale networks. SUMMARY

[0011] The first object of the present application is to provide a gene regulatory network topology identification method based on random variation Bayesian, and the present application identifies the topology structure of the gene regulatory network through a calculation method based on the maximization of the lower bound function J of Bayesian inference in random optimization.

[0012] The object of the present application is achieved by the technical scheme, and the specific steps are as follows:

[0013] 1) Data acquisition: acquiring a gene expression data set to be identified, the gene expression data set satisfying a Gaussian noise signal of normal distribution;

[0014] 2) Model construction: reconstructing the gene regulatory network model by adopting a dynamic structure function DSF based on a state space model;

[0015] 3) Model parameter estimation: estimating the parameters in the gene regulatory network model by adopting a method based on random variation Bayesian;

[0016] 4) Network topology structure identification: identifying the topology structure of the gene regulatory network by adopting a forward selection method, updating the network model by using an ARD variable, selecting the model structure by using the lower bound function J, and drawing a topology graph of the gene regulatory network.

[0017] Further, the specific structure of the reconstructed gene regulatory network model in step 2) is as follows:

[0018] Based on the state space model, the gene regulatory network is modeled by a dynamic structure function DSF, and the i-th gene can be described as:

[0019]

[0020] In the formula, y j (t) is the measured expression value of the i-th gene at t time, y j (t) is the measured expression value of the j-th gene at t time, Q ij is the regulatory effect of the j-th gene on the i-th gene, n is the number of genes, i, j [1, n], m is the number of inputs of the gene regulatory network, P ij is the regulatory effect of the j-th input on the i-th gene, u j (t) is the j-th input of the gene regulatory network, e i (t) is the measurement noise, Q ij and Pij is an element in Q and P, Q and P are transfer function matrices;

[0021] If the transfer function matrix Q is 0, it means that the ith gene is not regulated by the jth gene, and if the transfer function matrix P is 0, it means that the ith gene is not affected by the jth input, which is represented as no input in the gene regulatory network topology graph.

[0022] Further, the specific steps for estimating the parameters in the gene regulatory network model in step 3) are as follows:

[0023] 3-1) Convert the identification of the topology of the gene regulatory network into the estimation of the parameters of the gene regulatory network, denoted as q(Θ), which is equivalent to calculating the maximum value of the lower bound function J of the gene regulatory network parameters;

[0024] 3-2) Calculate the maximum of the lower bound function J based on the Bayesian variational method of stochastic optimization;

[0025] 3-3) Estimate the parameters of the gene regulatory network by the stochastic variational Bayesian method.

[0026] Further, the specific method for calculating the maximum value of the lower bound function J in step 3-1) is as follows:

[0027] The measured expression value of the ith gene can be described as:

[0028]

[0029] where h ij is the impulse response of Q ij in the model representation of the ith gene, l ij is the impulse response of P ij in the model representation of the ith gene, and Z is a sufficiently large positive integer used to make |h ij , Z|≈0, |l ij , Z|≈0;

[0030] Set the noise to follow a Gaussian distribution with a mean of zero, i.e.:

[0031] e i (t)~N(0,τ -1 ) (3)

[0032] where τ is the precision parameter of the Gaussian distribution of the measurement noise;

[0033] Set the prior distribution of τ to be a Gamma distribution, i.e.:

[0034]

[0035] where a0, b0 are two hyperparameters;

[0036] h ij and l ij The prior distribution is:

[0037]

[0038]

[0039] In the formula, and Represents the TC core, α j and β j Automatic correlation analysis of ARD variables to determine gene regulatory network topology;

[0040] Set α j and β j The prior distribution is a conjugate prior gamma distribution:

[0041] p(α j )=p(β j )=G(τ|a0,b0) (8)

[0042] kernel matrix and Each element in the array is defined as follows:

[0043]

[0044]

[0045] In the formula, γ 1，j and γ 2，j It is a constant of the kernel matrix;

[0046] The parameter is then expressed as ARD variables are represented as α = [α1, ..., α] n ], β=[β1,…,β n ], let Y i =[y i (1), ..., y i (N)] T ;

[0047] Network recognition based on Bayesian methods is the maximization of the likelihood function:

[0048]

[0049] In the formula, <·> represents the expectation operator, and q<·> represents the posterior distribution;

[0050] Let Θ=[h, l, τ, α, β] T Then the log-likelihood function of the output of the i-th gene in the gene regulatory network is:

[0051]

[0052] In the formula, Let KL divergence be the KL divergence. Let J be the lower limit of the likelihood function. Calculating q(Θ) is equivalent to calculating the maximum value of the lower limit J, i.e.:

[0053]

[0054] Furthermore, the specific method for maximizing the lower bound function J using the Bayesian variational method based on stochastic optimization in step 3-2) is as follows:

[0055] The lower bound function for ξ is:

[0056]

[0057] In the formula, ξ∈Θ, Θ -ξ ={Θ / ξ}, This represents the target expectation in terms of q(Θ), where const. indicates a constant.

[0058] The maximum value of J in terms of q(ξ) is given by:

[0059]

[0060] Its prior distribution p(ξ) and posterior distribution q(ξ) have conjugate properties and belong to the exponential family, that is:

[0061]

[0062] ln q(ξ)=ln s(ξ)+λ T t(ξ)-a g (λ) (17)

[0063] In the formula, λ is the natural parameter, t(ξ) is the sufficient statistic, and s(ξ) and a g (λ) represent the basic metric and the logarithmic regularizer, respectively;

[0064] The lower bound function J for q(ξ) can be rewritten as:

[0065]

[0066] In the formula:

[0067]

[0068]

[0069] The natural gradient of J(λ) has the following simple form:

[0070]

[0071] In the formula, the natural gradient is defined as the product of the standard gradient and the inverse of the Fisher information matrix of q(ξ), that is:

[0072]

[0073] The noise estimate of the natural gradient is obtained by subsampling the identification data, as follows:

[0074]

[0075] In the formula, t′ is uniformly sampled from {1, 2, ..., N}, and |t′| is the length of the sampled data;

[0076] If all gene nodes in a gene regulatory network are set independently, then:

[0077]

[0078] Based on the stochastic gradient algorithm, λ k It will converge to the optimum, which is updated by the following formula:

[0079]

[0080] λ k The convergence condition is:

[0081] ∑ k ρ k =∞, ∑ k (ρ k )2<∞ (26).

[0082] Furthermore, the specific method for estimating gene regulatory network parameters using the stochastic variational Bayesian method in step 3-3) is as follows:

[0083] 3-2-1) Initialize the posterior distribution q(Θ) and parameter γ 1，j and γ 2，j Learning rate ρ, parameter k;

[0084] 3-2-2) Update q(h) and q(l):

[0085] q(h ij The prior distribution of is a Gaussian distribution, and its exponential family model is:

[0086]

[0087] In the formula, λ1 is the natural parameter, and s1(h ij ) = 1, sufficient statistic t1(h ij ) is defined as Where Vec(·) is a vector consisting of the row vectors of a matrix;

[0088] The noise natural gradient of the lower bound function J(λ1) is:

[0089]

[0090] In the formula:

[0091]

[0092]

[0093] In the formula, h -j ={h\h ij}, φ j =[y i (t-1), ..., y i (tZ)] T , ψ j =[u j (t-1), …, u j (tZ)] T , Φ -j ={Φ\φ j};

[0094] The natural parameter λ1 is updated to:

[0095]

[0096] q(h) ij The distribution is Gaussian, with mean and covariance as follows:

[0097]

[0098]

[0099] In the formula, iVec is defined as the operator that converts a vector into a matrix. It is the inverse operation of Vec, and end represents the last index of the vector.

[0100] Similarly, q(l ij The exponential family model of ) is given by the following formula:

[0101]

[0102] The noise natural gradient of the lower bound function J(λ2) is:

[0103]

[0104] In the formula:

[0105]

[0106]

[0107] The natural parameter λ2 is updated to:

[0108]

[0109] Its mean and covariance are:

[0110]

[0111]

[0112] 3-2-3) Update q(α) j ) and q(β) j ):

[0113] q(α j The distribution is gamma, i.e.:

[0114]

[0115] In the formula, s3(α) j )=1,t3(α) j )=[α j ,lnα j ] T The natural gradient of the lower bound function J(λ3) is:

[0116]

[0117] q(α j Updated to:

[0118]

[0119] α j The mean is calculated as follows:

[0120]

[0121] Similarly, q(β) j It is also a gamma distribution, that is... The update rules are as follows:

[0122]

[0123] In the formula:

[0124]

[0125] 3-2-4) Update q(τ):

[0126] The posterior distribution of q(τ) is a gamma distribution, that is:

[0127]

[0128] The lower bound function of q(τ) is denoted by J(λ5), and its noise natural gradient is:

[0129]

[0130] q(τ) is updated to:

[0131]

[0132] 3-2-5) Update γ 1，j and γ 2，j :

[0133] In each iteration maximizing the lower bound, these two values ​​can be updated by solving two simple optimization problems, γ. 1，j The lower bound function is:

[0134]

[0135] In each iteration, γ 1，j and γ 2，j Updated to:

[0136]

[0137] 3-2-6) Calculate the lower bound J k If |J k -J k-1 If |<∈, where ∈ is the threshold, then stop; otherwise, jump to step 3-2-2).

[0138] 3-2-7) All parameters are estimated to their expected values.

[0139] Furthermore, the specific method for identifying the topology of the gene regulatory network using forward selection in step 4) is as follows:

[0140] The ARD variable is used to update the network model, and the lower bound function J is used to select the model structure, that is, the structure of the transfer matrices Q and P is determined by the lower bound function:

[0141] 4-1) Identify the structure of Q: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable α;

[0142] 4-1-1) ARD variable α j (j = 1, 2, ..., n) sorted in ascending order:

[0143]

[0144] 4-1-2) Let R = 1, then the gene regulatory network model is:

[0145]

[0146] 4-1-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R ;

[0147] 4-1-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrix Q; otherwise, R = R + 1 and jump to step 4-1-3.

[0148] 4-2) Identify the structure of P: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable β;

[0149] 4-2-1) ARD variable β j (j = 1, 2, ..., m) sorted in ascending order:

[0150]

[0151] 4-2-2) Let R = 1, then the gene regulatory network model is:

[0152]

[0153] 4-2-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R .

[0154] 4-2-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrices Q and P; otherwise, R = R + 1 and jump to step 4-2-3.

[0155] The second objective of this invention is to provide a gene regulatory network topology identification system based on stochastic variational Bayes, which includes:

[0156] The data acquisition unit is used to acquire the gene expression dataset to be identified, which is a Gaussian noise signal that follows a normal distribution.

[0157] The model building unit uses the dynamic structure function DSF based on the state space model to model and reconstruct the gene regulatory network model;

[0158] The model parameter estimation unit uses a stochastic variational Bayesian method to estimate the parameters in the gene regulatory network model.

[0159] The network topology identification unit uses a forward selection method to identify the topology of the gene regulatory network, updates the network model using ARD variables, and selects the model structure using a lower bound function J to draw the gene regulatory network topology diagram.

[0160] A third objective of the present invention is to provide a storage medium storing a plurality of instructions which are applicable to a processor for loading to execute any of the methods described above.

[0161] Because of the adoption of the above technical solution, the present invention has the following advantages:

[0162] 1. This application proposes to update the posterior distribution of global variables based on natural gradient noise from partial data, which greatly reduces the computational complexity.

[0163] 2. Compared with the full Bayesian inference method, this application can significantly reduce the amount of computation and is suitable for the identification of gene regulatory networks under large-scale data.

[0164] Other advantages, objectives, and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination, or may be learned from practice of the invention. The objectives and other advantages of the invention can be realized and obtained from the following description and claims. Attached Figure Description

[0165] The accompanying drawings of this invention are described below.

[0166] Figure 1 This is a schematic diagram of the method flow of the present invention.

[0167] Figure 2 This is a schematic diagram of the method flow of the present invention.

[0168] Figure 3 This is a topology diagram of the gene regulation network of the ring network type of the present invention;

[0169] Figure 4 This is a schematic diagram of gene regulatory networks identified from gene expression data. Detailed Implementation

[0170] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0171] like Figure 1 The method for identifying the topology of gene regulatory networks based on stochastic variational Bayes, as shown, includes the following specific steps:

[0172] 1) Data acquisition: Gene expression data of the gene regulatory network GRN to be identified were acquired from the database of the DREAM4 platform. The gene expression data is a Gaussian signal that meets the standard normal distribution.

[0173] 2) Model Construction: A dynamic structure function (DSF) based on a state-space model is used to model and represent the gene regulatory network in step 1); the specific method is as follows:

[0174] Based on the state-space model, the gene regulatory network is modeled using a dynamic structure function (DSF), and the i-th gene can be described as:

[0175]

[0176] In the formula, y i (t) represents the measured expression value of gene i at time t, y j (t) The measured expression value of gene j at time t, Q ij Let P represent the regulatory effect of gene j on gene i, where n is the number of genes, i, j ∈ [1, n], m is the number of inputs to the gene regulatory network, and P is the number of genes. ij For the regulatory effect of the j-th input on the i-th gene, u j (t) represents the j-th input to the gene regulatory network, e i (t) represents the noise level being measured, Q ij and P ij These are elements in Q and P, where Q and P are transfer function matrices.

[0177] 3) Model parameter estimation: The parameters in the model are estimated using a method based on stochastic variational Bayes. The specific steps are as follows:

[0178] 3-1) The identification of the gene regulatory network topology is transformed into the estimation of gene regulatory network parameters, denoted as q(Θ). The estimation of gene regulatory network parameters is equivalent to calculating the maximum value of its lower bound function J. The specific method is as follows:

[0179] The measured expression value of the i-th gene can be described as:

[0180]

[0181] In the formula, h ij For the i-th gene model, Q represents... ij The impulse response, l ij Let P be the representation of the i-th gene model. ij The impulse response, where Z is a sufficiently large positive integer, is used to make |h ij Z|≈0, |l ij Z|≈0;

[0182] Assume the noise follows a Gaussian distribution with a mean of zero, i.e.:

[0183] e i (t)~N(0,τ) -1 (55)

[0184] In the formula, τ is the accuracy parameter of the Gaussian distribution of the measured noise;

[0185] Let the prior distribution of τ be a gamma distribution, that is:

[0186]

[0187] In the formula, a0 and b0 are two hyperparameters;

[0188] h ij and l ij The prior distribution is:

[0189]

[0190]

[0191] In the formula, and Represents the TC core, α j and β j Automatic correlation analysis of ARD variables to determine gene regulatory network topology;

[0192] Set α j and β j The prior distribution is a conjugate prior gamma distribution:

[0193] p(α j )=p(β j )=G(τ|a0,b0) (59)

[0194] kernel matrix and Each element in the array is defined as follows:

[0195]

[0196]

[0197] In the formula, γ 1，j and γ 2，j It is a constant of the kernel matrix;

[0198] The parameter is then expressed as ARD variables are represented as α = [α1, ..., α] n ], β=[β1,…,β n ], let Y i =[yi (1), ..., y i (N)] T ;

[0199] Network recognition based on Bayesian methods is the maximization of the likelihood function:

[0200]

[0201] In the formula, <·> represents the expectation operator, and q<·> represents the posterior distribution;

[0202] Let Θ=[h, l, τ, α, β] T Then the log-likelihood function of the output of the i-th gene in the gene regulatory network is:

[0203]

[0204]

[0205] In the formula, Let KL divergence be the KL divergence. Let J be the lower limit of the likelihood function. Calculating q(Θ) is equivalent to calculating the maximum value of the lower limit J, i.e.:

[0206]

[0207] 3-2) The lower bound function J is maximized using a Bayesian variational method based on stochastic optimization. The specific method is as follows:

[0208] The lower bound function for ξ is:

[0209]

[0210] In the formula, ξ∈Θ, Θ -ξ ={Θ / ξ}, This represents the target expectation in terms of q(Θ), where const. indicates a constant.

[0211] The maximum value of J in terms of q(ξ) is given by:

[0212]

[0213] Its prior distribution p(ξ) and posterior distribution q(ξ) have conjugate properties and belong to the exponential family, that is:

[0214]

[0215] ln q(ξ)=ln s(ξ)+λ T t(ξ)-a g (λ) (68)

[0216] In the formula, λ is the natural parameter, t(ξ) is the sufficient statistic, and s(ω and a) is the s-statistic. g (λ) represent the basic metric and the logarithmic regularizer, respectively;

[0217] The lower bound function J for q(ξ) can be rewritten as:

[0218]

[0219] In the formula:

[0220]

[0221]

[0222] The natural gradient of J(λ) has the following simple form:

[0223]

[0224] In the formula, the natural gradient is defined as the product of the standard gradient and the inverse of the Fisher information matrix of q(ξ), that is:

[0225]

[0226] The noise estimate of the natural gradient is obtained by subsampling the identification data, as follows:

[0227]

[0228] In the formula, t′ is uniformly sampled from {1, 2, ..., N}, and |t′| is the length of the sampled data;

[0229] If all gene nodes in a gene regulatory network are set independently, then:

[0230]

[0231] Based on the stochastic gradient algorithm, λ k It will converge to the optimum, which is updated by the following formula:

[0232]

[0233] λ k The convergence condition is:

[0234] ∑ k ρ k =∞, ∑ k (ρ k )2<∞ (77).

[0235] 3-3) Estimate the parameters of the gene regulatory network using the stochastic variational Bayesian method. The specific steps are as follows:

[0236] 3-3-1) Initialize the posterior distribution q(Θ) and parameter γ 1，j and γ 2，j Learning rate ρ, parameter k;

[0237] 3-3-2) Update q(h) and q(l):

[0238] q(h ij The prior distribution of is a Gaussian distribution, and its exponential family model is:

[0239]

[0240] In the formula, λ1 is the natural parameter, and s1(h ij ) = 1, sufficient statistic t1(h ij ) is defined as Where Vec(·) is a vector consisting of the row vectors of a matrix;

[0241] The noise natural gradient of the lower bound function J(λ1) is:

[0242]

[0243] In the formula:

[0244]

[0245]

[0246] In the formula, h -j ={h\h ij}, φ j =[y i (t-1), ..., y i (tZ)] T , ψ j =[u j (t-1), …, u j (tZ)] T , Φ -j ={Φ\φ j};

[0247] The natural parameter λ1 is updated to:

[0248]

[0249] q(h) ij The distribution is Gaussian, with mean and covariance as follows:

[0250]

[0251]

[0252] In the formula, iVec is defined as the operator that converts a vector into a matrix. It is the inverse operation of Vec, and end represents the last index of the vector.

[0253] Similarly, q(l ij The exponential family model of ) is given by the following formula:

[0254]

[0255] The noise natural gradient of the lower bound function J(λ2) is:

[0256]

[0257] In the formula:

[0258]

[0259]

[0260] The natural parameter λ2 is updated to:

[0261]

[0262] Its mean and covariance are:

[0263]

[0264]

[0265] 3-3-3) Update q(α) j ) and q(β) j ):

[0266] q(α j The distribution is gamma, i.e.:

[0267]

[0268] In the formula, s3(α) j )=1,t3(α) j )=[α j ,lnα j ] T The natural gradient of the lower bound function J(λ3) is:

[0269]

[0270] q(α j Updated to:

[0271]

[0272] α jThe mean is calculated as follows:

[0273]

[0274] Similarly, q(β) j It is also a gamma distribution, that is... The update rules are as follows:

[0275]

[0276] In the formula:

[0277]

[0278] 3-3-4) Update q(τ):

[0279] The posterior distribution of q(τ) is a gamma distribution, that is:

[0280]

[0281] The lower bound function of q(τ) is denoted by J(λ5), and its noise natural gradient is:

[0282]

[0283] q(τ) is updated to:

[0284]

[0285] 3-3-5) Update γ 1，j and γ 2，j :

[0286] In each iteration maximizing the lower bound, these two values ​​can be updated by solving two simple optimization problems, γ. 1，j The lower bound function is:

[0287]

[0288] In each iteration, γ 1，j and γ 2，j Updated to:

[0289]

[0290] 3-3-6) Calculate the lower bound J k If |J k -J k-1 If |<∈, where ∈ is the threshold, then stop; otherwise, jump to step 3-2-2).

[0291] 3-3-7) All parameters are estimated to their expected values.

[0292] 4) Network topology identification: The forward selection method is used to identify the topology of the gene regulatory network. The ARD variable is used to update the network model, and the lower bound function J is used to select the model structure. The specific method is as follows:

[0293] The specific method for performing network structure identification using the forward selection method in step 4) is as follows:

[0294] The ARD variable is used to update the network model, and the lower bound function J is used to select the model structure, that is, the structure of the transfer matrices Q and P is determined by the lower bound function:

[0295] 4-1) Identify the structure of Q: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable α;

[0296] 4-1-1) ARD variable α j (j = 1, 2, ..., n) sorted in ascending order:

[0297]

[0298] 4-1-2) Let R = 1, then the gene regulatory network model is:

[0299]

[0300] 4-1-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R ;

[0301] 4-1-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrices Q and P; otherwise, R = R + 1 and jump to step 4-1-3.

[0302] 4-2) Identify the structure of P: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable β;

[0303] 4-2-1) ARD variable β j (j = 1, 2, ..., m) sorted in ascending order:

[0304]

[0305] 4-2-2) Let R = 1, then the gene regulatory network model is:

[0306]

[0307] 4-2-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R .

[0308] 4-2-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrices Q and P; otherwise, R = R + 1 and jump to step 4-2-3.

[0309] In this invention, the basic idea for determining the topology is to use ARD variables to update the network model and use J to select the model structure; whereby ARD variables can be used to evaluate the strength between nodes: if α j Or β j Close to the initial value, which means that in this case, h ij or l ij If the variance matrix is ​​close to zero, and Q or P is close to 0, then it can be inferred that the i-th node is not controlled by the j-th node or the input.

[0310] A forward selection method is employed. First, a fully connected network model is estimated using stochastic variational inference from step 3). The resulting ARD variables (α and β) are ordered in ascending order. Nodes are added to the model sequentially. In each step, a specific network model is used for stochastic VI estimation until the likelihood function no longer increases.

[0311] 5) Data simulation: The typical gene regulatory network is simulated using the methods described above. The specific methods are as follows:

[0312] Set the hyperparameter of the gamma distribution to a0 = 10. -2 and a0 = 10 -4 An adaptive strategy was used to design a stochastically optimized learning rate, which was set as follows:

[0313]

[0314] Where ρ0, k s ε∈[0.5,1] are the setting parameters. Accuracy (ACC), True Positive Rate (TPR), and Precision (PREC) were used to evaluate the performance of the recognition algorithm. All simulation experiments were run on a Macmini desktop with a 3.2GHz CPU.

[0315] 5-1) Identify gene regulatory networks of random network type, the specific methods are as follows:

[0316] Consider a gene regulatory network with 5 nodes and 5 input terminals. The state-space model of this gene network is as follows:

[0317] x(t+1)=Ax(t)+Bu(t)+e(t) (106)

[0318] y(t)=Cx(t) (107)

[0319] In the formula, A and B are two sparse stable matrices, and the non-zero elements of A and B are given by the following formula:

[0320] A(3,1)=0.0096, A(2,2)=0.4403, A(3,2)=0.1853, (108)

[0321] A(4,4)=-0.0756, A(5,4)=1.4980, B(2,1)=-1.3433, (109)

[0322] B(3,1)=-0.8883, B(4,1)=0.8221, B(3,2)=0.3222, (110)

[0323] B(5,2)=-1.1761, B(1,3)=0.8395, B(4,4)=0.4461, B(3,5)=-0.2217 (111)

[0324] The input is gene expression data represented by a Gaussian signal with a standard normal distribution. The process noise is also a Gaussian signal with the same variance as the input. Each experiment generates 1000 data points. To demonstrate the proposed recognition method based on random VI, the recognition of the first node is studied without loss of generality. Based on the proposed method, we only need to sample a subset of data points to update the model parameters. Let η represent the ratio of the number of samples to the number of data points in each iteration. If η = 1, the proposed method becomes a traditional VI. The lower bound curves for different η and different adaptive learning rates are shown below. Figure 2 As shown.

[0325] The random VI-based method converges slightly faster on all data (η=1) than on partial data (η≠1). However, the convergence results on partial data are similar to those of the full VI method; with a well-chosen learning rate, the convergence results on partial data are also similar to those of the full VI method. The advantage of the random VI method is its lower computational cost compared to the traditional VI method. Table 1 lists the computational costs for different sampling rates. Clearly, computational complexity increases with increasing sampling rate. Based on these results, we can see that by selecting appropriate sampling and learning rates, the proposed algorithm achieves good convergence and low computational cost with low computational complexity. This is crucial for identification using large-scale data.

[0326] Table 1 Calculation costs under different η values

[0327] η = 0.01 η = 0.05 η = 0.1 η = 0.2 η = 0.5 η = 1 149.82s 159.77s 172.19s 197.41s 276.06s 420.38s

[0328] To further test the proposed method, we generated two stable random sparse gene regulation networks (sparse density set to 0.3), with n = m = 5 and n = m = 10, respectively. The input and process noise were Gaussian signals with the same variance, fixed at 1, and 1000 data points were collected for each experiment. Fifty randomized experiments were conducted, and the average performance index is shown in Table 2. Full Bayesian inference has significant advantages in network identification and is the state-of-the-art method for network identification. Therefore, we compared the proposed method with the full Bayesian method. As shown in Table 2, random VIs ​​with sampling rates of 5% and 10% can achieve similar performance to traditional inference on network topology, but with significantly reduced time.

[0329] Table 2 Performance of Random VI on Different Networks

[0330]

[0331]

[0332] 5-2) Identify circular network-type gene regulatory networks, using the following specific methods:

[0333] Consider two extremely sparse gene regulatory networks with a large number of nodes. Figure 3 The structure of this network is shown in Figure 5. Each node is controlled by the next node and the next input (the last node is controlled by the first node and the next input). Given the non-zero elements of A: A(i, i+1) = 1 ≤ i ≤ n-1, A(n, 1). The input matrix B is the same as A. The input signal and process noise are the same as in Figure 5-1.

[0334] The random VI method was tested on large-scale networks with n=m=50 and n=m=100. 10,000 sampling points were generated for each network. The sampling rate for random VI was chosen to be 0.05. The results are listed in Table 3. The results show that all network nodes and the input matrix were correctly identified. This clearly demonstrates that the proposed method is effective for these extremely sparse networks. Furthermore, with η=0.05, obtaining a network with 50 nodes took nearly 100 hours, while obtaining a network with 100 nodes took approximately 560 hours. This method was not compared with the fully Bayesian method because the computational complexity of the fully Bayesian method is too high to be tested on our platform. However, simulation results indicate that due to the lower computational cost of the proposed method, it is possible to apply Bayesian methods to such large-scale gene regulatory networks.

[0335] Table 3 Performance indices of two circular gene regulatory networks

[0336] n = m = 50 n = m = 100 PREC 1 1 TPR 1 1 ACC 1 1 Time (h) 99.58 561.47

[0337] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0338] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0339] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0340] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0341] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for identifying the topology of gene regulatory networks based on stochastic variational Bayes, characterized in that, The specific steps are as follows: 1) Data acquisition: Acquire the gene expression dataset to be identified, which is a Gaussian noise signal that follows a normal distribution; 2) Model Construction: The gene regulatory network model is reconstructed using dynamic structure function DSF based on the state-space model; 3) Model parameter estimation: The parameters in the gene regulatory network model are estimated using a method based on stochastic variational Bayes. 4) Network topology identification: The forward selection method is used to identify the topology of the gene regulatory network, the ARD variable is used to update the network model, and the lower bound function J is used to select the model structure and draw the gene regulatory network topology diagram. The specific structure of the reconstructed gene regulatory network model in step 2) is as follows: Based on the state-space model, the gene regulatory network is modeled using a dynamic structure function (DSF), and the i-th gene is described as follows: In the formula, y i (t) represents the measured expression value of gene i at time t, y j (t) represents the measured expression value of the j-th gene at time t, Q ij Let P represent the regulatory effect of gene j on gene i, where n is the number of genes, i, j ∈ [1, n], m is the number of inputs to the gene regulatory network, and P is the number of genes. ij For the regulatory effect of the j-th input on the i-th gene, u j (t) represents the j-th input to the gene regulatory network, e i (t) represents the noise level being measured, Q ij and P ij These are elements in Q and P, where Q and P are transfer function matrices; If the transfer function matrix Q is 0, it means that the i-th gene is not regulated by the j-th gene. If the transfer function matrix P is 0, it means that the i-th gene is not affected by the j-th input. In the gene regulation network topology diagram, this is represented by the node having no input. Step 3) involves estimating the parameters in the gene regulatory network model. The specific steps are as follows: 3-1) The identification of the topology of the gene regulatory network is converted into the estimation of the parameters of the gene regulatory network, denoted as q(Θ). The estimation of the parameters of the gene regulatory network is equivalent to calculating the maximum value of its lower bound function J. 3-2) Calculate the lower bound function J maximization using a Bayesian variational method based on stochastic optimization; 3-3) Estimate the parameters of the gene regulatory network using the stochastic variational Bayesian method; The specific method for identifying the topology of the gene regulatory network using forward selection in step 4) is as follows: The ARD variable is used to update the network model, and the lower bound function J is used to select the model structure, that is, the structure of the transfer matrices Q and P is determined by the lower bound function: 4-1) Identify the structure of Q: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable α; 4-1-1) ARD variable α j (j = 1, 2, ..., n) sorted in ascending order: 4-1-2) Let R = 1, then the gene regulatory network model is: 4-1-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R ; 4-1-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrix Q; otherwise, R = R + 1 and jump to step 4-1-3. 4-2) Identify the structure of P: For fully connected gene regulatory network models, use the stochastic variational Bayesian method in step 3) to estimate the ARD variable β; 4-2-1) ARD variable β j (j = 1, 2, ..., m) sorted in ascending order: 4-2-2) Let R = 1, then the gene regulatory network model is: 4-2-3) For the gene regulatory network model, the lower bound is calculated using the stochastic variational Bayesian method in step 3), denoted as J. R ; 4-2-4) If J R <J R-1 If the condition is met, stop and output the transfer function matrices Q and P; otherwise, R = R + 1 and jump to step 4-2-3.

2. The method for identifying the topology of gene regulatory networks based on stochastic variational Bayes as described in claim 1, characterized in that, The specific method for calculating the maximum value of its lower bound function J in step 3-1) is as follows: The measured expression value of the i-th gene is described as follows: In the formula, h ij For the i-th gene model, Q represents... ij The impulse response, l ij Let P be the representation of the i-th gene model. ij The impulse response, where Z is a sufficiently large positive integer, is used to make |h ij，Z |≈0,|l ij，Z |≈0; Assume the noise follows a Gaussian distribution with a mean of zero, i.e.: e i (t)~N(0,τ -1 ) (3) In the formula, τ is the accuracy parameter of the Gaussian distribution of the measured noise; Let the prior distribution of τ be a gamma distribution, that is: In the formula, a0 and b0 are two hyperparameters; h ij and l ij The prior distribution is: In the formula, and Represents the TC core, α j and β j Automatic correlation analysis of ARD variables to determine gene regulatory network topology; Set α j and β j The prior distribution is a conjugate prior gamma distribution: p(a j )=p(β j )=G(τ|a0,b0) (8) kernel matrix and Each element in the array is defined as follows: In the formula, γ 1，j and γ 2，j It is a constant of the kernel matrix; The parameter is then expressed as ARD variables are represented as α = [α1, ..., α] n ], β=[β1,…,β n ], let Y i =[y i (1), ..., y i (N)] T ; Network recognition based on Bayesian methods is the maximization of the likelihood function: In the formula, <·> represents the expectation operator, and q<·> represents the posterior distribution; Let Θ=[h, l, τ, α, β] T Then the log-likelihood function of the output of the i-th gene in the gene regulatory network is: In the formula, Let KL divergence be the KL divergence. Let J be the lower limit of the likelihood function. Calculating q(Θ) is equivalent to calculating the maximum value of the lower limit J, i.e.:

3. The method for identifying the topology of gene regulatory networks based on stochastic variational Bayes as described in claim 1, characterized in that, The specific method for maximizing the lower bound function J using the Bayesian variational method based on stochastic optimization in step 3-2) is as follows: The lower bound function for ξ is: In the formula, ξ∈Θ, Θ -ξ ={Θ / ξ}, Let J represent the target expectation in terms of q(ξ), where const. indicates a constant; and let J have the maximum value in terms of q(ξ): Its prior distribution p(ξ) and posterior distribution q(ξ) have conjugate properties and belong to the exponential family, that is: lnq(ξ)=lns(ξ)+λ T t(ξ)-a g (l) (17) In the formula, λ is the natural parameter, t(ξ) is the sufficient statistic, and s(ξ) and a g (λ) represent the basic metric and the logarithmic regularizer, respectively; the lower bound function J for q(ξ) is rewritten as: In the formula: The natural gradient of J(λ) has the following simple form: In the formula, the natural gradient is defined as the product of the standard gradient and the inverse of the Fisher information matrix of q(ξ), that is: The noise estimate of the natural gradient is obtained by subsampling the identification data, as follows: In the formula, t′ is uniformly sampled from {1,2,…,N}, and |t′| is the length of the sampled data; If all gene nodes in a gene regulatory network are set independently, then: Based on the stochastic gradient algorithm, λ k It will converge to the optimum, which is updated by the following formula: λ k The convergence condition is: ∑ k r k =∞,∑ k (r k ) 2 <∞ (26).

4. The method for identifying the topology of gene regulatory networks based on stochastic variational Bayes as described in claim 1, characterized in that, The specific method for estimating gene regulatory network parameters using the stochastic variational Bayesian method in step 3-3) is as follows: 3-2-1) Initialize the posterior distribution q(Θ) and parameter γ 1,j and γ 2,j Learning rate ρ, parameter k; 3-2-2) Update q(h) and q(l): q(h ij The prior distribution of is a Gaussian distribution, and its exponential family model is: In the formula, λ1 is the natural parameter, and s1(h ij ) = 1, sufficient statistic t1(h ij ) is defined as Where Vec(·) is a vector consisting of the row vectors of a matrix; The noise natural gradient of the lower bound function J(λ1) is: In the formula: where h -j ={h \ h ij}, φ j =[y i (t - 1),..., y i (t - Z)] T , ψ j =[u j (t - 1),..., u j (t - Z)] T , Φ -j ={φ \ φ j}; The natural parameter λ1 is updated to: q(h ij The distribution is Gaussian, with mean and covariance as follows: In the formula, iVec is defined as the operator that converts a vector into a matrix. It is the inverse operation of Vec, and end represents the last index of the vector. Similarly, q(l ij The exponential family model of ) is given by the following formula: The noise natural gradient of the lower bound function J(λ2) is: In the formula: The natural parameter λ2 is updated to: Its mean and covariance are: 3-2-3) Update q(α) j ) and q(β) j ): q(α j The distribution is gamma, i.e.: In the formula, s3(α) j )=1,t3(α) j )=[α j ,lnα j ] T The natural gradient of the lower bound function J(λ3) is: q(α j Updated to: α j The mean is calculated as follows: Similarly, q(β) j It is also a gamma distribution, that is... The update rules are as follows: In the formula: 3-2-4) Update q(τ): The posterior distribution of q(τ) is a gamma distribution, i.e.: The lower bound function of q(τ) is denoted by J(λ5), and its noise natural gradient is: q(τ) is updated to: 3-2-5) Update γ 1,j and γ 2,j : In each iteration maximizing the lower bound, these two values ​​are updated by solving two simple optimization problems, γ. 1,j The lower bound function is: In each iteration, γ 1,j and γ 2,j Updated to: 3-2-6) Calculate the lower bound J k If |J k -J k-1 If |<∈, where ∈ is the threshold, then stop; otherwise, jump to step 3-2-2). 3-2-7) All parameters are estimated to their expected values.

5. A storage medium, characterized in that, The storage medium stores a plurality of instructions which are adapted for loading by a processor to execute the method of any one of claims 1 to 4.

Images