Approximate inference of gene regulatory network models from RNA-Seq time series data

Background Inference of gene regulatory network structures from RNA-Seq data is challenging due to the nature of the data, as measurements take the form of counts of reads mapped to a given gene. Here we present a model for RNA-Seq time series data that applies a negative binomial distribution for the observations, and uses sparse regression with a horseshoe prior to learn a dynamic Bayesian network of interactions between genes. We use a variational inference scheme to learn approximate posterior distributions for the model parameters. Results The methodology is benchmarked on synthetic data designed to replicate the distribution of real world RNA-Seq data. We compare our method to other sparse regression approaches and find improved performance in learning directed networks. We demonstrate an application of our method to a publicly available human neuronal stem cell differentiation RNA-Seq time series data set to infer the underlying network structure. Conclusions Our method is able to improve performance on synthetic data by explicitly modelling the statistical distribution of the data when learning networks from RNA-Seq time series. Applying approximate inference techniques we can learn network structures quickly with only moderate computing resources.


Background
Methods for the inference of gene regulatory networks from RNA-Seq data are currently not as mature as those developed for microarray datasets. Normalised microarray data posses the desirable property of being approximately normally distributed so that they are readily amenable to various forms of inference, and in the literature many graphical modelling schemes have been explored that exploit the normality of the data [1][2][3][4][5][6][7][8][9].
The data generated by RNA-Seq studies on the other hand present a more challenging inference problem, as the data are no longer approximately normally distributed, and before normalisation take the form of non-negative integers. In the detection of differential expression in RNA-Seq data, negative binomial distributions have been applied [10][11][12][13], providing a good fit to the overdispersion typically seen in the data relative to a Poisson distribution. Following similar graphical modelling approaches as applied in the analysis of microarray data, it is natural to consider Poisson and negative binomially distributed graphical models. Unfortunately in many cases Correspondence: t.thorne@reading.ac.uk Department of Computer Science, University of Reading, Reading, UK when applying graphical modelling approaches with Poisson distributed observations, only models that represent negative conditional dependencies are available, or inference is significantly complicated due to lack of conjugacy between distributions [14]. Poisson graphical models have been applied successfully in the analysis of miRNA regulatory interactions [15,16], but we might expect to improve on these by modelling the overdispersion seen in typical RNA-Seq data sets with a negative binomial model.
One specific case of interest in the analysis of RNA-Seq data is the study of time series in a manner that takes into account the temporal relationships between data points. Previous work in the literature has developed sophisticated models for the inference of networks from microarray time series data [4,5], but whilst methods have been developed for the analysis of differential behaviour in RNA-Seq time series [17,18], little attention has been given to the task of learning networks from such data. Although existing nonparametric methods applicable to time series may be applied [19,20], these were not specifically designed for application to RNA-Seq data, and also require time consuming approaches such as Markov Chain Monte Carlo schemes. There are also existing information theoretic methods, for example those of [21], but again these were designed for application to microarray data, and are not designed for time series data and learning of directed networks.
Here we present a method for the inference of networks from RNA-Seq time series data through the application of a Dynamic Bayesian Network (DBN) approach, that models the RNA-Seq count data as being negative binomially distributed conditional on the expression levels of a set of predictors. Whilst there has been work applying negative binomial regularised regression approaches in the literature [22], here we specifically consider the problem of learning networks from RNA-Seq data, and apply the horseshoe prior [23,24], that has been shown to have advantages in robustness and adaptivity over other regularisation methods.

Dynamic Bayesian Networks
In a DBN framework [25], considering only edges between time points, we can model a sequence of observations using a first order Markov model, where the value of a variable at time t is dependent only on the values of a set of parent variables at time t − 1. This is illustrated in Fig. 1 and can be written as where Pa(i) is the set of parents of variable i in the network. In our case we wish to model the expression level of a gene conditional on a set of parent genes that have some influence on it. To learn the set of parent variables of a given gene, it is possible to perform variable selection in a Markov Chain Monte Carlo framework, proposing to add or remove genes to the parent set in a Metropolis-Hastings sampler. Another option is to consider all possible sets of parent genes as suggested in [20]. However for even modestly sized sets of genes (e.g. 50) this can be computationally expensive, and so instead we consider applying a sparse regression approach to learn a set of parents for each gene. This approach considers the contribution of all possible parent genes in a regression framework but encourages sparsity in the coefficients so that only a small set are non-zero.

Sparse negative binomial regression
Given data D consisting of M columns and L rows, with columns corresponding to genes and rows to time points, we seek to learn a parent set for each gene. To do so we can employ a regularised regression approach that enforces sparsity of the regression coefficients, and only take predictors (genes) whose coefficients are significantly larger than zero as parents. To simplify the presentation, below we consider the regression of the counts for a single gene i, y = D 2:L,i , conditional on the counts of the remaining The matrix X is supplemented with a column vector 1 to include a constant term in the regression. Where there are multiple replicates for each time point these can be adjusted appropriately.
The counts y t are then modelled as following a negative binomial distribution with mean exp (Xβ) t and dispersion ω, where β is a vector of regression coefficients β w and a constant term β c . The model for a gene i is then where we have applied the NB2 formulation of the negative binomial distribution, and s t is a scaling factor for each sample to account for sequencing depth. The s t can be estimated from the data by considering the sum of counts for each sample, or by the more robust approach of [11] where the median of ratios is used. We place a straightforward normal prior on β c and to enforce sparsity of the β w we apply a horseshoe prior [23,24], assuming that β w ∼ N (0, ζ 2 w ), and placing a half-Cauchy prior on the ζ 2 w , Then as in [24] we set a prior on τ that allows the degree of shrinkage to be learnt from the data An example of the sparsity induced in the β w can be seen in figure 8 in Appendix 2. Finally we place a gamma prior on the dispersion parameter ω. This gives a joint probability (subsuming the model parameters into θ) of

Variational Inference
We now apply a variational inference [26][27][28][29][30] scheme to learn approximate posterior distributions over the model parameters. In a Bayesian setting variational inference aims to approximate the posterior p(θ|x) with a distrubtion q(θ). To do so we attempt to minimise the Kullback-Leibler (KL) divergence between the two, defined as As the KL divergence is bounded below by zero, it follows from 9 that and so we can define a lower bound on the logarithm of the model evidence as To make the problem of minimising the KL divergence tractable we can consider a mean field approximation where the posterior is approximated by a series of independent distributions q(θ i ) on some partition of the parameters, Under the mean field assumption it can be shown that to minimise the KL divergence between p(θ|x) and q(θ), or equivalently to maximise the model evidence lower bound (or ELBO) L(q), the optimal form for each q(θ i ) is where the expectation is over the remaining q(θ j =i ). In many cases this formalism is sufficient to derive a coordiante ascent algorithm to maximise the ELBO where the variational parameters of theq(θ i ) are updated iteratively. Unfortunately in our model the optimal distribution q for the regression coefficients β w does not have a tractable solution. However following [31] we can sidestep this problem by applying non-conjugate variational message passing [32], and we can then derive approximate posterior distributions for each of the model parameters following a straightforward parameter update scheme. The full set of variational updates are given in Appendix 1.
Considering our model as a graphical model as in Fig. 2, we can decompose the terms of E q j =i log p(θ, x) in Eq. 14 into those dependent on θ i by considering the neighbours of θ i . Then we can rewrite Eq. 14 as where Ch i denotes the children of node i in the graphical model. Considering each term on the right hand side of Eq. 15 as a message from another variable in the graphical model it is possible to deriveq in the conjugate exponential family as in [33]. In the non-conjugate case, the messages can be approximated as in [32], derived for the negative binomial model in [31].

Synthetic data
We apply our method to the task of inferring directed networks from simulated gene expression time series. The time series were generated by utilising the GeneNetWeaver [34] software to first generate subnetworks representative of the structure of the Synthetic count data were generated by constructing a negative binomial DBN model as in Eq. 2 corresponding to the generated subnetworks with randomised parameters β sampled from a mixture of equally weighted N (0.3, 0.1)  and N (−0.3, 0.1) distributions. The initial conditions and mean and dispersion parameters were randomly sampled from the empirically estimated means and dispersions of each gene from a publicly available RNA-seq count data set from the recount2 database [35] (accession ERP003613) consisting of 171 samples from 25 tissues in humans [36]. This was done so as to simulate the observed distributions of RNA-Seq counts in a real world data set.
We compare our approach against the Lasso as implemented in the lars R package [37], and the Gaussian regularised regression method in the glmnet R package [38]. For these methods the count data was first normalised, either by transforming the counts by the empirical cumulative distribution function of the data and subsequently mapping these to the quantiles of a N (0, 1) distribution, or by applying the rlog function of the DESeq2 R package [13] to normalise the counts. We also applied the regularised Poisson regression method implemented of the glmnet R package to the count data, and the mpath R package [22] that performs penalised negative binomial regression. Finally we also applied a multinomial regularised regression from the glmnet R package to discretised data that were binned into 4 distinct levels by quantiles, to give a discrete DBN model. The degree of regularisation was in each case selected using cross validation as implemented in the respective software packages.
In Figs. 3 and 4 we show the partial area under the receiver operating curve (AUC-ROC) with a cutoff of 0.95 and corrected to fit the range 0 to 1, and area under the precision recall curve (AUC-PR), as calculated by the PRROC R package [39], and Matthews Correlation Coefficient (MCC), for the various methods to be benchmarked. For the MCC, edges were predicted as those where zero Fig. 3 Boxplots of partial AUC-ROC, AUC-PR, and MCC for our method (Nb) and the competing methods benchmarked when learning directed networks of 25 nodes from synthetic data, for 5 subnetworks sampled from the S. cerevisiae gene regulatory network Fig. 4 Boxplots of partial AUC-ROC, AUC-PR, and MCC for our method (Nb) and the methods benchmarked when learning directed networks of 50 nodes from synthetic data, for 5 subnetworks sampled from the S. cerevisiae gene regulatory network was not contained in the 95% credible interval of the corresponding regression coefficients, and for the Lasso and glmnet methods, non-zero coefficients were taken as predicted edges. As the count data were generated by a stochastic model, we repeated benchmarking on each network 5 times with resampled negative binomial means and dispersions and simulated count data. Running time for our algorithm was under 10 minutes for the 50 node networks considered.
For networks of 25 nodes in Fig. 3, our method shows an improved performance over the competing methods in terms of the AUC-PR, and also in terms of the MCC. Although the distinction between the approaches is less marked for the AUC-ROC, this is to be expected as the simulated biological network structures have far fewer (< 10%) true positives than true negatives, a situation in which the AUC-ROC does not distinguish performance as well as AUC-PR [39,40].
As can be seen in Fig. 4 performance for larger networks of 50 nodes is also improved over competing methods in terms of AUC-PR and MCC. For the competing methods, quantile normalisation for the Lasso and glmnet appear to outperform normalisation using the rlog function of DESeq2. As the only other method applying a negative binomial distribution, mpath is the closest method to our approach, but it appears that the application of the horseshoe shrinkage prior delivers improved performance. It is clear that, as might be expected, taking into account the distributional properties observed in RNA-Seq data improves on the performance of methods based on assumptions that do not hold for RNA-Seq count data.

Neural progenitor cell differentiation
To illustrate an application of our model to a real world RNA-Seq data set, we consider a publicly available RNA-Seq time course data set available from the recount2 database [35], accession SRP041159. The data consist of RNA-Seq counts from neuronal stem cells for 3 replicates over 7 time points after the induction of neuronal differentiation [41]. To select a subset of genes to analyse we performed a differential expression test between time points using the DESeq2 R package [13], and selected the 25 genes with the largest median fold-change between time points that were also differentially expressed between all time points.
Applying our method and selecting edges with a posterior probability > 0.95 produced the network shown in Fig. 5, where it can be seen that there are four genes (MCUR1, PARP12, COL17A1, CDON) acting as hubs, suggesting these genes may be important in neuronal differentiation. Within the network MCUR1 appears to influence the transcription of a large number of genes with many outgoing edges, whilst PARTP12, COL17A1 and CDON have both incoming and outgoing edges. This may suggest a more fundamental role of MCUR1 in controlling neuronal differentiation.
For each node we also calculate the betweenness centrality, which is the fraction of the total number of shortest paths between nodes in the network that pass through a given node. This gives a measure of the importance of a node in the network, as nodes with a larger betweeness centrality will disrupt more paths within the network if deleted, and act as bottlenecks that connect modules within the network. Looking at the betweenness centrality of each node it appears that PARP12 and CDON, and to a lesser extent COL17A1, are important carriers of information within the network. Of these genes playing a central role in the network, CDON has been shown to be promote neuronal differentiation through the activation of p38MAPK pathway [42,43] and inhibition of Wnt signalling [44], whilst MCUR1 is known to bind to MCU [45], that in turn has been shown to influence neuronal differentiation [46].

Discussion and conclusions
We have developed a fast and robust methodology for the inference of gene regulatory networks from RNA-Seq data that specifically models the observed count data as being negative binomially distributed. Our approach outperforms other sparse regression based methods in learning directed networks from time series data.
Another approach to network inference from RNA-Seq data could be to further develop mutual information based methodologies with this specific problem in mind. Mutual information based methods have the benefit of being independent of any specific model of the distribution of the data, and so could help sidestep issues in parametric modelling of RNA-Seq data. However this comes at the cost of abandoning the simplifying assumptions that are made by applying a parametric model that provides a reasonable fit to the data, and presents challenges of its own in the reliable estimation of the mutual information between random variables.

Appendix 1: Variational inference
From the results in [31] the model can be written as a Poisson-Gamma mixture, so that and the horseshoe prior on β represented using a mixture of inverse gamma distributions,

Mean field approximation
The mean field approximation of the posterior is then The variational updates for λ t can be derived as logq(λ t ) = E q log p(y t |λ t )p(λ t |X t , β, ω) + const.
and due to the properties of the log-normal distribution where is the covariance matrix of β underq. As derived in [31], applying non-conjugate variational message passing,q(β) ∼ N (μ, ) and the variational update for β is and for the dispersion ω we apply numerical integration as described in [31]. Then for the horseshoe prior on β, the variational updates are logq ζ 2 w = E q log p β w |ζ 2 w p ζ 2 w + const.
Appendix 2: Supplemental Figures  Fig. 6 Metrics calculated for networks of 25 nodes separated by individual network structure for the 5 different networks considered. Each bar plot corresponds to 5 simulated data sets from a single network structure