Robust inference of the context specific structure and temporal dynamics of gene regulatory network

Regulatory networks are specific to exogenous context

Since gene or protein in the network can acquire different functions under different exogenous conditions or phenotypic contexts, the network constituents and interactions would be context dependent, i.e., regulatory network is composed of different genes and proteins under different conditions and their interactions are specifically determined by the context. For instance, the oncogenes such as MYC, ERK, and WNT are known to down-regulated in normal cells. Up-regulation of oncogenes leads to a series of modification in signal transduction and transcription regulation, which interfere with the controlled cell death and eventually turn the normal cell into cancer cell. Therefore, the revelation of context-specific differences of regulatory network in cancer is essential for understanding the process and mechanism of cancer. In fact, the context-dependent regulation accounts for the heterogeneity of diverse cell phenotypes.

Regulatory networks are time varying according to the endogenous cell states

While the long run behaviour of regulatory networks is often stationary, the transient characteristic is intrinsically time varying for they need to respond and adapt themselves to the new exogenous condition. Temporal dynamics of regulatory networks is highly prevalent to many important processes such as cell cycle, cell differentiation, apoptosis, and viral infection to account for different endogenous cellular states. Even under a fixed exogenous condition, the regulation of genes and protein are carried out in a temporal, sequential manner. The cell cycle is known to have multiple stages and transitions in between and are triggered and guided by temporal changes of regulatory networks. Some of the transcription factors such as Fos are known to be early genes because they regulate commonly at the early stage of a biological event. It is also well recognized that the regulatory pathway does not always persist over all times. Such temporal cascade of regulation is essential for elucidating the mechanism leading to different phenotypes. The two key characteristics of time-varying networks include 1) the temporal topology of regulatory networks, i.e., the network constituents and wiring could change with time, and 2) time-varying regulatory impact. It is important and necessary to understand both aspects of temporal regulation when reconstructing the time varying networks.

Existing network models and approaches are insufficient for inferring context-specific time-varying networks.

Reconstruction of regulatory networks based on high throughput data is currently a research focus in systems biology. In the past, a large number of different network architectures have been proposed. One of the most popular architectures is the molecular-based network architecture, where nodes are individual genes or proteins and links denote gene-gene, gene-protein, or protein-gene regulations. This architecture aims to reveal the detailed genomic interaction and is capable of uncovering transcription regulation, signalling, and metabolic pathways. However, the quality of the inferred molecular-based network is greatly restricted by the limitation of data availability and data quality; in practice, precise inference of large scale networks at a molecular level is very difficult and not robust due to the large number of the unknowns and limited and noisy genomics or proteomics data. The reconstructed networks normally have extremely high false positive links, making the subsequent analysis and interpretation at a molecular level difficult and unreliable. A common remedy to this problem is to take a global view of the networks and study the features such as scale free and modularity. Although this approach has been successful in revealing global properties of molecular networks, it is still unable to identify pathways as intended by this architecture. An alternative and equally popular architecture is the module networks [2]. Since nodes in a module network are clusters of genes, or modules, the inference task for module networks is less demanding and higher quality inference can be achieved. However, module networks cannot directly uncover functional pathways even with high quality data and thus are often difficult to interpret biologically. Post-processing based on gene ontology (GO) or pathway enrichment is required to reveal the possible functions of each modules. However, the module networks cannot reveal beyond the prediction that a set of functions of a modular would regulate another set in another module. This prediction is useful for pointing out the possible further studies to follow; it is nevertheless far from sufficient to uncover important pathways and transcription regulations. Apparently, a new architecture that is robust to data noise and capable of revealing directly pathways will be much more appealing.

Given a network structure, a mathematical or statistical model needs to be imposed to model the dynamics of regulation. Many models have been proposed for this purpose including ordinary differential equations [3], (probabilistic) Boolean networks [4], [5], information theory based models [6]. However, the principled assumption underlying these models is overwhelmingly stationary; that is the dynamic models do not change over time. While the long run behaviour of regulatory networks is reasonably considered stationary, the transient characteristic is intrinsically time varying for they need to respond and adapt themselves to meet the need of varying conditions. As a result, the inferred networks based on stationary assumption conceal the intricate temporal regulatory patterns. Despite the prevalent temporal nature of regulatory, limited work exists to model and infer the time varying network and they include regime state space model [7], temporal hidden Markov model [8], the wavelet dynamic vector autoregressive method [9], and context-sensitive probabilistic Boolean networks [10]. These models all involve a high degree of parametric complexity that hinders robust reconstruction of large scale networks. For instance, context-sensitive probabilistic Boolean networks in [10] and regime state space model [7] involved only 10 and 47 genes, respectively. Alternatively, module networks were assumed in [7–9] to increase the inference robustness but nevertheless reduce the interpretability of result. Some remarkable results about systems level reconstruction of both static and transient transcription regulatory networks of yeast for different cell states were reported in [11] .The networks inference utilized a large collection of prior knowledge and a simple model of differential expression in microarray. However, specific time regimes for cell states of interest need to be given; the approach cannot further uncover by itself endogenous temporal patterns within the given regime. In addition, the adopted gene networks architecture is sensitive to data noise and prone to producing false positive links.

Aside from modelling network dynamics, investigation of the context-specific networks is given increasing attention, although the main body of existing work still concerns the reconstruction of context-independent networks. To this end, [12] explicitly introduced a context node in the Bayesian network model to model different context condition. GO annotations were applied in [13] to the protein-protein interaction (PPI) network of yeast to dissect the context (function)-specific sub¬networks. Again, in [11], large scale transcription networks under different exogenous conditions were constructed for yeast. Despite some success in identifying context-specific network motifs and transcription factors, the results were not robust and unreliable at the molecular scale due to the molecular-based network architecture..

To sum up, making robust statistical inference from diverse types of genomic data to discover context-specific and temporal pathways of biological function presents a formidable challenge. The current network architectures and approaches are insufficient for this goal. The main issue is that the network architecture, model, and inference approach cannot produce robust network inference that can be easily understood.

Towards a pathway-centric network model and functional enrichment inference approach

To overcome the problems of existing network architecture, model, and inference approaches, we advocate a new pathway-centric network (PathRNet) and an enrichment inference approach. The nodes of the PathRNet are regulators such as transcription factors (TFs) ,pathways including signaling pathways and metabolic pathways, and the links bear specific functional meaning, i.e., a link going from TF to pathway defines the transcription regulation of pathway by TF and a link from signaling pathway to TF suggests signal transduction induced modification of TF. The advantages of this model architecture are in four folds. First, this network design is geared directly to biologists, who normally focus on specific pathways and regulators when designing experiments and seeking answers to phenotypic outcomes. Secondly, PathRNet speaks the language of biologists and no interpretation is even needed for the results. Unlike in most of gene or protein networks, where links only imply a degree of association and hence further interpretable is needed, in PathRNet, the linkage such as Signaling pathway A -> TF B -| Metabolic pathway C reads naturally as Signaling pathway A activates TF B, which further inhibits metabolic pathway C. Thirdly, past research has produced considerable data about pathways and transcription regulation, and it is prudent and necessary to understand their regulatory roles in different context before new pathways and regulations can be inferred. Fourthly, PathRNet is more robust to data noise. Since nodes are pathways instead of individual genes, PathRNet inherits the robustness of module networks; given these advantages, we believe that PathRNet can adequately address the aforementioned deficiency of the existing network architectures.

To robustly infer the context specific structure and its temporal dynamics of PathRNet, we propose the Pathway Enrichment Analysis of Temporal Regulatory Networks or PATTERN. The key of PATTERN is the method of the enrichment based time dependent clustering algorithm, ECTDISA [14], which identifies temporal clusters of pathways and regulated gene sets that are enriched in a context-specific time series microarray data. Once again, the inference in PATTERN moves away from the gene-centric approach to the pathway or gene-set enrichment based concept. The gene set enrichment (GSE) concept formulated by the GSEA algorithm [15] has gained great success in differential analysis of gene expression data and also extended its applications to drug effect prediction and clustering of coregulated gene. The main advantages of GSE approach are 1) the inference is robust to excessive noise on one or a few genes in the gene set, and 2) the result has straightforward biologically implications since gene sets are mostly pathways or coregulated genes. Recalling the previously mentioned problems of current inference approaches for network reconstruction, we see that these advantages of GSE are the exact answers to these problems. However, GSE has not been used for network construction. It is thus our objective to develop an enrichment approach for network reconstruction. Coupled with PathRNet, the GSE approach will provide added robustness to inference, enabling reliable discovery of network temporal dynamics from noisy microarray data of commonly small sample size.

Overview of our approach

The block diagram of our approach PATTERN is shown in Fig. 1. It has a modular structure including four modules. In module 1, pathway data and TF regulated gene sets data from the existing knowledge database are assembled to construct a generic PathRNet. To construct the network, only two types of directed links are allowed, i.e., a link from a signaling pathway node to a TF node is placed if TF is in the downstream of the signaling pathway and a link from a TF node to a signaling or a metabolic pathway is placed if the TF regulated gene set contains genes in the pathway. The resulting PathRNet is a generic network that represents the prior knowledge of regulation under variety of different contexts. In Module 2, context-specific microarray data is introduced and the context-specific PathRNet is obtained by a functional enrichment strategy. Particularly, ECTDISA [14] is applied to determine the enrichment of each node in the generic PathRNet in the time series microarray data and the nodes that are not enriched will be deactivated. Furthermore, ECTDISA also determines the temporal expression pattern (the expression trend and time span) of each node. In Module 3, regulatory PathRNet is constructed to depict the regulatory relationship of nodes based on the expression pattern obtained in Module 2. Inference is carried out to determine the presence of links and their regulatory impact. As a result, a (context-specific) regulatory PathRNet will be obtained. In Module 4, temporal regulatory impact of network during each time period is further analyzed, resulting in a temporal landscape of PathRNet dynamics.

We discuss next the proposed methods for implementing each module.

Generic PathRNet Construction

A generic PathRNet is first constructed based on the existing knowledge database. The prior knowledge databases are collected from NCI-NC [16] ,Transfac [17] and MSigDB [15]. A generic PathRNet is constructed with each node corresponding to a signaling pathway, a transcription factor or a metabolic pathway; and an edge is placed:

a) From a signaling pathway to TFs if TFs are the downstream of the signaling pathways regulate transcription factors.

b) From a TF to signaling or metabolic pathways if TF regulated at least one genes in the pathway.

In the end, a prior network is obtained, which illustrates a generic picture of the interactions between pathways and transcription factors inside a human cell.

Noted that in our approach, there are additional steps to refine the constructed generic PathRNet by removing redundant nodes and edges:

In the context-specific PathRNet construction module, all the nonfunctional nodes and their edges were removed, leaving only the enriched functional pathways /TFs and edges between them.

In the regulatory PathRNet construction module, Bayesian information criterion (BIC) was applied for model selection, and all the edges that are not chosen were removed from the constructed regulatory PathRNet.

In the dynamic PathRNet construction module, (dynamic) edges that indicate no or nominal contributions were further removed from the constructed temporal PathRNet.

Since all the following steps remove redundant nodes/edges, we hope keep as many associations as possible when constructing the generic PathRNet to ensure a high sensitivity of reconstruction. Thus, an edge is placed between a transcription factor (TF) and a pathway as long as the TF regulates at least one gene in the Pathway (the pathway and the TF target gene sets share at least one gene). However, the degree of overlapping between the pathway genes and the TF target gene set might itself indicates the strength of their association. Although this information is not used in our algorithm, it is possible to include this information in further extension to better estimate the regulatory impacts of transcription factors.

Context specific PathRNet construction

To capture the functional components of generic PathRNet and construct a context-specific PathRNet, expression profile is then incorporated. The goal is to identify coregulated pathways enriched by the gene expression profile and their associated time period. To this end, we apply the ECTDISA [14] to the microarray gene expression data. ECTDISA can not only cluster the coregulated genes that are functionally most enriched by a functional annotation database such as pathway databases but also retrieve the timing of the co-regulating clusters.

Particularly, when applying the ECTDISA, a measurement needs to be specified to gauge the coregulation of genes. Since Pearson's correlation is not robust towards noise, and Euclidean distance is not suitable to measure correlation, we proposed a modified Euclidean distance measurement ρ to penalize smaller fold changes without losing the robustness to noise, which is defined as follows:

where x, y are the expression profile of a functional node (pathway or TF), and α is used to penalize smaller fold changes.

Since small fold changes are more vulnerable to noise, the parameter α was chosen to penalize smaller fold change to a reasonable amount. It should not be too small or too large. An experiment is developed on real data to show the sensitivity of the parameter, and result is evaluated in terms of A score and C score, which are defined in [18] for performance evaluation purpose. In general, A score indicates whether all the identified modules are corresponding to true modules; C score indicates whether all the true modules are correctly identified by clustering algorithm. Since our goal is to choose a proper α to uncover all the modules from the data, we are mainly interested in acquiring the highest C score; but at the same time, A score should not be too low. The result is shown in Fig. 2.

It can be seen from the figures that C score reaches its peak when α = 0.3 and at the same time A score is relatively high. This is the value we choose when the approach was applied to real data. To achieve the best result, the use of different α are recommended.

An example of clusters produced by ECTDISA is shown in Fig. 3. A number of the most significantly enriched signaling pathways, transcription factors and metabolic pathways are selected as the functional components in the prior network, and the context specific PathRNet is constructed by retaining only the functional components.

Regulatory PathRNet construction

Note that for now the edges in the functional networks are still defined by the prior knowledge. With the expression data, the edges can be further refined to reflect the regulatory status including up- or down- regulation and its weight. To this end, the expression profiles of the nodes are first determined. For the real case, since even the genes in the same pathways can behave differently, i.e., some are up regulated while some others are down regulated, it possible that the same pathway are enriched in both the up-regulated module and the down-regulated module. To capture the main trend of pathway, we consider the most significantly enriched trend in a pathway by using ECTDISA algorithm. Specifically, for a signaling pathway or a metabolic pathway (nodes labelled with "S-" and "P-" in the network), the expression profile is estimated by averaging all the pathway genes that are most enriched in their temporal trends. Since the trends of their expression profiles are very similar, the averaging expression of those genes is a good estimate of the main trend of the pathway. For example, if at 1 hour, the expressions of 5 genes that are involved in the same pathways are 2, 2.1, 2.2, 0, -3; and the first 3 gene are enriched in a module identified by ECTDISA, then the fold change of the corresponding pathway at 1 hour will be . Similarly, the expression fold change at other time points can be calculated. Note that, since the expression is calculated by averaging along the gene dimension, the expression fold change of different genes are acquired from exactly the same experimental condition; moreover, the microarray data has been properly normalized [19]

If a pathway is enriched in a temporal transcription module, then the pathway is only active during the inferred time period that the temporal module is active; otherwise its expression is 0. For a transcription factor, the expression fold change profile is inferred by averaging the expression fold changes of its encoding genes.

Next, the regulatory impact is inferred based on the expression of nodes. Let y denote the expression profiles of a node and X = [x₁, x₂,...,x _k ] the expression profiles of its parent nodes. Then, a linear regulatory model is introduced as:

y = X · β + ε

where β represents regulatory weights and ε is white Gaussian noise. The goal is then to determine the functional parent nodes and the corresponding regulatory weights β. The problem is a variable selection problem and BIC is adopted for the solution. Under the assumption that the model error is normally distributed, the formula for BIC is:

where n is the number of data points in x _i , k is the number of parameters to be estimated, and the residual sum of squares (RSS) is the sum of squares of residuals, which defines the discrepancy between the data and our estimation model:

RSS=|| y - ŷ||²

where, is an estimate of y calculated based on the least squares estimate:

Temporal PathRNet construction

In the last module, the dynamic regulatory network is constructed from the basal network by considering the dynamic contribution of the regulatory nodes at each sample time points. Specifically, the dynamic change of the expression level of a node at a sample time t can be expressed as:

where, x _t = [x_{1, t} , x_{2, t} ,…,x _p,t ], and P is the number of parent nodes of y. The dynamic contribution of all the parent nodes is (x _t – x _{t –1})∙β, and the dynamic contribution of a parent node from t – 1 to t is: . Then, the dynamic regulatory networks can be obtained by refining the edges at each t. Particularly, if is small, the corresponding edge can be removed indicating no or nominal contribution. Otherwise, the width and the color of the edge will be defined by to illustrate the regulatory influence of the node onto a child node during a particular sample time period. Although small regulatory contribution does lead to edge removal and thus change of network structure, it is mainly the regulatory time period inferred by the ECTDISA that defines the dynamics of the network structure.

Generic PathRNet

The generic PathRNet, as shown in Fig. 9, was constructed by integrating the existing knowledge derived from NCI-NC, Transfac and MSigDB [15–17], which includes 324 signaling pathways, 945 transcription factors, and 259 metabolic pathways. The constructed generic PathRNet consists of 281 nodes representing TFs, signaling / metabolic pathways and 1354 linkages, illustrating a generic picture of all the known interactions between pathways and transcription factors for different human cells under various conditions.

Context-specific PathRNet revealed by ECTDISA

A time series microarray data of KSHV infection of HUVEC [19] was introduced and the context-specific PathRNet (Fig. 10) for KSHV infection was reconstructed. The raw data consist of a total of 22283 features (probe set IDs) sampled at time t = [0,1,3,6,10,16,24,36,54,78] (hour) after infection. An intensity filter (the intensity of a gene should be above 100 in at least 1 sample) and a variance filter (the inter-quartile range of log2-intensities should be at least 0.2) were then applied to select 7159 differentially expressed features, which are further merged into 5650 Entrez gene database UIDs by taking the maximum value of all corresponding probe set IDs.

ECTDISA[18] was then applied to the time series data to obtain the coregulated gene clusters that are most enriched by the nodes of the generic PathRNet and at the same time retrieve the timing of coregulation.

Regarding how many nodes should be kept in the context-specific PathRNet, a node is considered to be significantly enriched in a module when the enrichment log p-value is larger than a predefined threshold. The threshold is chosen to be dependent on the noise level of the data. In general, a stricter threshold may result in too few nodes, i.e., low recall, and thus misses many true network connections; a relatively loose threshold may result in too many nodes, i.e., high false positive, and includes many false connections. As experiment was developed to show the sensitivity of the threshold and the result is shown in Fig. 11.

Specifically, PATTERN was applied to a simulated dataset which consists of 400 genes 10 samples. The PathRNet embedded is shown in Fig. 4. Different p-value thresholds are applied and compared. It can be seen from Fig. 11 that when the threshold is too strict (-40) or too loose (-2.5) the proposed algorithm cannot retrieve the true network structure. The result is better when a reasonable threshold is chosen. To this end, we used a p-value threshold of 0.01 to keep relatively more nodes in the network.

By keeping the only enriched or functional nodes (table 1) of the generic PathRNet, the context-specific PathRNet for KSHV infection was obtained. Since a number of enriched nodes are independent, those nodes are removed from the network. The reconstructed PathRNet consists of 38 nodes and 208 regulations (Fig. 10). Each node corresponds to a signaling pathway (S-), transcription factor (TF-) or a metabolic pathway (P-). The node enriched in the same cluster is highlighted by the same color and also indicated by the class ID "C#". Each directed edge represents the regulatory relationship.

Regulatory and temporal PathRNet

The temporal PathRNets for the KSHV infection of HUVEC were constructed to reveal the dynamics in regulation at each sampling time of the microarray experiment.

To better illustrate the result, a sub-network of related to p38 MAPK signaling pathway (S7) is acquired by keeping the nodes that are within the 2 nodes away from S-7. 9 frames of this sub-network are shown in Fig. 12. Please refer to Fig. 3 for the expression profiles of the involved clusters.

First of all, the time varying nature of the network structure and regulatory impact are revealed in Fig. 13, 14, 15, 16, 17, 18, 19, 20, 21. The color of the nodes indicates the expression fold changes during the time period specified in Fig. 13, 14, 15, 16, 17, 18, 19, 20, 21 (Red represents up-regulation, and blue represents down-regulation). A directed linkage is determined by β _p in Fig. 12, which is the regulation coefficient, and by in Fig 13, 14, 15, 16, 17, 18, 19, 20, 21, which are the dynamic contributions of the parent nodes to the child node during the time period specified. (Red represents up-regulation, and green represents down-regulation, Please refer to Additional File 1, Appendix for the annotation of each node.)

In the Fig 13, 14, 15, 16, 17, 18, 19, 20, 21, the gray color of the nodes and links indicates that they are not functional. As shown by these figures, the colors of nodes and links vary over time, revealing the time-dependent characteristics of the network. Regarding the variation of the network structure, two aspects of the algorithm are important. First, the time-dependent module defines the start and the end of an enriched connection in the module. For instance, since P-128 is enriched in a temporal transcription module that is not defined between 0 to 3 h post infection (hpi), P-128 does not function in the first 2 frames of the time varying network (Fig. 13 and 14) and thus the nodes and links are in gray. However, they become active after 3 hpi (Fig. 15-16 and 20-21). Secondly, when the regulatory contribution of a parent node is sufficiently small, the parent node and the link are also marked gray to indicate nominal regulation. For example, while TF-75 regulates P-123, -126, -129 and -100 at different time points; it only regulates some or even none at some instances, say, between 3-6 hpi, or it only regulates P-129 and P-126 at other time points. Furthermore, it is easy to observe that the contributions for many nodes vary at different time points. Interestingly, a positive feedback loop exists between TF-12 and S-7.

Next, interesting regulatory cascades can be revealed by incorporating the annotations with the time varying networks. In Fig. 13 from 0-1hpi, the p38 MAPK pathway (S7) is up-regulated, which also induces the TF STAT1 (TF-61). This suggests that in the early stage of infection KSHV activates p38 MAPK pathway. Indeed, the activation of the p38 MAPK pathway is shown to be essential for KSHV infection and inhibition of the p38 MAPK pathway reduced KSHV infectivity [21]. Induction of STAT1 likely indicates the activation of cellular innate immune response to viral infection [22]. At 3 hpi in Fig. 14, it is suggested that STAT1 also activates SMAD2/3 signaling (S47), which is in agreement with the existing knowledge that activation MAPK pathway enhances Smad2 transcriptional activity [22]. In addition, TF CREB1 (TF-75) is also activated during both 1 and 3 hpi, which is expected since it is at the downstream of the MAPK pathway. Interestingly, the entire cascade from p38 MAPK to CREB1 becomes downregulated at 6 hpi. Following the first wave of activation of the cascade, we have also observed a second wave of activation of the p38 MAPK pathway and TF STAT1 at 10 hpi. The signal is again transduced to SMAD2/3 at 16 hpi accompanying the decay of the p38 MAPK pathway and TF STAT1, before being further transduced to CREB1 at 24 hpi. Activation of this second wave of the p38 MAPK pathway has indeed been described in [21, 23].