Reverse-engineering of TF network and TFBS information

A number of techniques are available to uncover the topology of the TF network—the networks of complex reactions and interactions in the cell that control transcript levels [12]. One strategy is to use principals of reverse-engineering and use gene expression data to infer regulatory interactions [13, 14]. Various reverse-engineering methods can reduce the dimensionality of the classic combinatorial search problem and utilize genome sequence data to enhance the sensitivity and specificity of predictions. However, they have difficulties in describing regulatory control by mechanisms other than TFs. Reverse-engineering of TF networks in the lower eukaryotes has been well developed [15–17]. However, the problems in mapping the regulatory mechanisms in cells of higher eukaryotes have made such global studies either impossible or impractical. Some recent studies have begun to address this issue [18–20], but have tended to focus only understanding which TFs bind to which genes—not looking in detail at the nature of the TF/TFBS interaction. A recent study [21] identified key biological features in transcriptional changes, however this method has difficulties in inferring the dynamics of the interactions. Furthermore, TF concentrations were not considered during the identification of the features.

To date, various reverse-engineering methods can reduce the dimensionality of the reverse-engineering problem and utilize genome sequence data to enhance the sensitivity and specificity of predicted interactions. However, they have difficulties in describing regulatory control by mechanisms other than TFs. To address this issue TFBSs information is required to complement the gene expression data. We used a list of 132,654 TFBSs between 20,920 genes and 174 TFs that had been identified by searching an alignment of five mammal species for conserved 5’ and 3’ regions [22]. Connectivity data is notorious for high false positive rates; however, our connectivity data is robust against the problem because it extracts binding information from well conserved upstream regions. A more detailed explanation is addressed in the Methods and Results section and a schematic diagram of the connectivity data is presented in Fig. 1.

Identifying regulation type by combining TFA and TFC analysis

Transcription factor activities (TFAs) are the intensity of the interactions between a certain transcription factor (TF) and its targets at a certain experimental point [23]. Thus, the estimated strength of TFAs between each TF and its target gene are useful to know which TF is acting on which gene at a given time point or experiment condition. However, simply knowing the regulatory activities under a single experimental condition provides limited information about the transcriptional network. To understand the mechanism of regulatory interactions, we developed a method that identifies statistically significant differences in TFAs under two different conditions. The significant differences indicate the changing level of TFAs between two conditions, so varying trends of TFAs in whole experimental process are easily detected and can be used to identify TF-specific regulatory patterns (up and down-regulation).

A highly concentrated TF induces more gene expressions rather than a lower concentrated TF. High-affinity binding sites induce the gene expressions at any level of TF concentration (TFCs), but low-affinity binding sites require high level of TF concentration for induction [24]. Thus, we might assume that TF concentration level is an important factor for investigating TFAs, and TFA investigation with considering TFC provides more reliable and accurate results closer to the complex reality of biology. To address this problem we proposed a probabilistic variational inference method to infer the concentration of each TF protein (TFC) and the regulatory intensities (TFA) of each TF and gene pair [4].

Aside from the method, there have been some notable attempts to infer TFAs based on integrating gene expression data and TFBSs information. The approaches use various well-known statistical inference techniques such as network component analysis [25], support vector machine [26], multivariate regression plus backward variable selection [27] and partial least squares [28]. However, the TFAs, which are inferred by these methods, do not contain any information on the strength and the sign of the physical interaction between a TF and its target genes. Moreover, the regulatory interactions can change easily in response to changing experimental conditions and over time. Since the methods are not fully probabilistic, they are not ideal for investigating the stochastic interactions. A linear regression based probabilistic method to model the full probability distribution of each TFA on each gene was developed [23]. The limitation of this method, however, is that it does not infer the TFAs and TFCs separately. This is a serious problem in subsequent analysis and prediction.

Transcription regulatory circuits and mathematical model

Typically, we have knowledge of \( \underset{\bar{\mkern6mu}}{\boldsymbol{e}} \) (from gene expression profiling experiments, such as microarray or RNA-seq) and would like to infer the set of TFC \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \), and TFA \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) giving rise to this signal. Given \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \) and \( \underset{\bar{\mkern6mu}}{\boldsymbol{e}} \) Sanguinetti et al. [4] then show how it is possible to solve for \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \) and \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) using a discrete time state space model (Eq. 1) with expectation-maximization (EM) algorithm. In the model, elements of the \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \) matrix indicate the concentration level of a given TF protein (TFC) at a specific time. Elements of the \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) matrix represent the regulatory intensity (TFA) between a given TF protein and its binding affinity to its target genes. The baseline expression level is the mean vector. The measurement noise \( \underset{\bar{\mkern6mu}}{\boldsymbol{v}} \) follows zero-mean i.i.d. Gaussian noise. To estimate the \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \) and \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) matrices, the model used posterior estimation of Bayer’s theorem. During this estimation, EM algorithm allowed the model to efficiently approximate the log likelihood. However, it is rare to have a complete knowledge of \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \) —we simply do not know the binding sites for all TFs in a typical higher eukaryotic cell. Recent experimental techniques, such as ChIP-chip and ChIP-seq can provide useful data to help construct the connection topology \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \) [31], however they have clear limitations if we are looking for a complete topology [32, 33]. A number of theoretical techniques are also available to uncover the connection topology [34–36]. The techniques generally use principals of reverse-engineering and use gene expression and genome sequence data to infer regulatory interactions.

Gene expression data

Gene expression datasets were downloaded from Gene Expression Omnibus (accession number GSE7342 for p38α and GSE36890 for STAT5) [37, 38]. The expression profiling data of GSE7342 dataset was normalized by the robust microarray average (RMA) method. The read counts of GSE36890 dataset was normalized to the reads per kilo-base of exon per mega-base of library size (RPKM).

Generating \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \), the connection topology

In this paper we have taken a conservative strategy for generating \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \) which looks at upstream region of genes that are well-conserved in multiple mammalian genomes. We used a published catalogue of common regulatory motifs that were overrepresented in gene upstream regions [22]. These motifs were identified by constructing genome-wide alignments for four mammalian species in promoter regions and 3’ UTRs relating to well-annotated genes from the RefSeq database. The same TFs were assumed to bind the same TFBSs in mice since the TFBS had been discovered in an alignment of human, mouse, rat and dog promoter regions. TFBS upstream of human 13,330 RefSeq genes were predicted. Mouse genes corresponding to the published list of human genes [22] were identified using Ensembl mouse gene annotation.

Estimation of statistically significant changes

The value of the Cutoff was chosen such that all differences at the 95% confidence interval were considered significant (±2 standard deviations). ±2SD limit is widely chosen as a normal limit because it fits well into two important categories: (1) confident interval and (2) testing hypothesis.

Gene Ontology (GO) analysis

GO analysis was performed by using DAVID [39]. The sets of genes showing significant changes identified in Eq. 2 were submitted to DAVID using the default parameters in order to obtain the GO term classifications of each gene. Our computational pipeline utilized the results to investigate the functionality of genes and their regulatory TFs. The detailed methods and the result figures of GO analysis are supplied in Additional file 1: Supplementary Text and Figure S1–S3.

Estimating the responses to perturbation of transcription networks

We have developed a strategy which used forward-engineering to construct the connection topology (see Fig. 1, Methods, and Additional file 1: Supplementary Text and Figure S1–S3), based on a previous study of regions upstream of genes conserved in multiple mammalian genomes [22]. The structure of this network of transcriptional regulatory interactions between TFs and the genes whose transcription they control is described by a binary matrix \( \boldsymbol{T}\ \in\ {\mathbf{\Re}}^{\boldsymbol{n}\times \boldsymbol{m}} \), where n is the number of TFs and m is the number of genes; An element (i, j) of the matrix is ‘1’ if TF i binds to the upstream control region of gene j, ‘0’ otherwise. We have then employed a mathematical model to integrate the connection topology data and a gene expression dataset from a higher eukaryote in which we are interested in modeling the changes that occur in the TF network in response to a change in the cellular environment (Fig. 2). Our approach could be seen as complementary to ‘Integrative methods’, as defined in [40], as it provides a strategy for creating an approximate connection topology if more detailed information is not available. The connection topology that is being used for this analysis contains many approximations and is certainly incomplete. However, it should be noted that we are looking at the differences between the models, for example between a wild-type and knock-out state, and those differences will be in parts of the model for which we do have data.

The results of our approach provide a set of TFs and their target genes which are related by significant up- or down-regulation in transcription. It provides a clear indication of the changes in TFA and TFC of TFs that are controlling transcriptional regulatory mechanisms in response to a specific stimulus. We therefore showed an “integrated” approach for network inference, based on a forward-engineered connection topology, can produce plausible and testable hypotheses about the responses to perturbation of transcription networks in higher eukaryotes.

Illustrating interpretable images of complex data

The visualization tools then make patterns apparent that would be difficult to detect in numerical data (Fig. 2). To distinguish regulation patterns between different experimental conditions, recognizing at a glance is important. However, computing results are formed in large numerical matrix, thus it is not only difficult to navigate through the whole matrix but also impossible to present the results in one page.

Figure 3 shows a graphical representation of the significant changes in TFA matrices \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) (n by m) and TFC vectors \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \) (n) obtained from this analysis. The patterns of the responses to perturbation of TF networks are readily observed in this single-shot image that presents approximately 2000 significant changes of varying TF activities on the 132,654 TFBSs after deleting p38α. In upper part of the plots, the TFs place in the functional group order. The genes, which have at least one significant interaction with TFs, locate in bottom part of the plots. A line in the plots presents a regulatory interaction (normalized TFA by TFC) between a TF and its target gene, and line color indicates a significant difference between the strengths of the regulatory interaction of two conditions. For example, we can easily find in visualized format (Fig. 3) that TF group three has distinct patterns (down-regulation at E13.5, up-regulation at E15.5) between two time points.

Modelling the changes of transcription factor network in p38α deficient mice

The computational pipeline as highlighted in Fig. 2 was applied to a published study of the effect p38α knock-out in mouse embryos [37]. This study developed four gene expression profiling datasets (Gene Expression Omnibus, accession number GSE7342) comprising of two time points at days 13.5 and 15.5 of embryonic development (E13.5 and E15.5) for p38α knock-outs and their wild-type controls. This data set was chosen for this study as it includes experimental measurements of gene expression in the wild-type and knock-out mice and showed that p38α deficient mice have significantly different phenotype. Thus, the experimental datasets were used as positive controls for our theoretical study.

From the TFA weight \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{W}}} \) and connection topology \( \underset{\bar{\mkern6mu}}{\underset{\bar{\mkern6mu}}{\boldsymbol{T}}} \) matrices, the average strength of TFA, S _f , for each TF in the datasets was calculated:

By comparing the average strengths between wild-type and knock-out mice, it is possible to see which of the TFs have significantly changed as a consequence of the removal of p38α.

Figure 4a, b show the changes in TFA strengths between wild-type and knock-out mice at E13.5 and E15.5. It can be seen that a number of TFs show a significant signal (>2 s.d.) in this data. These are shown with more detail in Table 1. Figure 4c, d show the inferred TFCs \( \underset{\bar{\mkern6mu}}{\boldsymbol{c}} \) obtained for the E13.5 and E15.5 time points. Again, from this graph it is possible to see that a number of TFs appear to be responding to the p38α status. These are shown in more detail in Table 2.

Trans name	TF group	Cis name	Regul-ation type	GO analy-sis	Known biological functions	Ref.
OLF1	1.Dev	EBF1	Down	Imm	NF-κB → EBF1 → B-cell development	[52]
GATA-X	1. Dev	GATA1 ~6	Down	-	p38 → HSp27 → GATA1 → Differentiation → GATA1 → IL9 → Asthma	[53, 54]
HNF3	1. Dev	FOXM1	Down	-	p38 → FOXM1 → Apoptosis	[45]
RSRFC4	5. Res	MEF2	Down	Apo, Dev Imm	p38 → MEF2 → Development → Differentiation	[46, 55]
NF-κB	6. Lat	NF-κB	Down	Imm	p38 → NF-κB → IFNγ → STAT1 → Development → ZEB1 → Immune	[56–58]
SOX9	7. Unk	SOX9	Down	Apop, Dev	p38 → SOX9 → Apoptosis	[59–62]
PITX2	1. Dev	PITX2	Up	Apo	p38 → PITX2 → Development → Apoptosis	[63, 64]
HFH4	1. Dev	FOXJ1	Up	-	FOXJ1 ─┤ NF-κB	[65]
FXR	3. Ste	FXR	Up	Apo, Dev	p38 → FXR → Apoptosis	[66, 67]
AP1	5. Res	c-JUN	Up	-	p38 ─┤JNK → c-JUN → Proliferation → Apoptosis	[37]
MYC	5. Res	MYC	Up	Apo, Dev	p38 ─┤MYC → Apoptosis	[68]

The TFs predicted to have significantly different behaviors between the wild-type and knock-out mice. “Trans Name”—the official gene symbol of the TF. “Cis Name”—the name given to the binding site of the TF. TFs were characterized into different functional groupings (see Fig. 3c for details: Dev—cell-type specific developmental TFs, Res—signal dependent resident nuclear factors, Lat—signal dependent latent cytoplasmic factors, Ste—signal dependent steroid receptor group, Unk—unknown). “Regulation Type”, the way in which the TF regulates its target genes in the absence of p38α. “GO Analysis”, provides more functional classification for the identified TFs (see Methods and Additional file 1: Figure S1 to S3 for details). The abbreviations of the GO terms are: Apo, Apoptosis; Dev, Developmental Process; Imm, Immune System Development. “Known Biological Functions” summarizes the findings from the recent biological literature as shown in the “References”

The boldface ones are the main node in Fig. 5

Trans name	TF group	Cis name	Level	GO analy-sis	Known biological functions	Ref.
AREB6	1.Dev	ZEB1	Down	-	p38 → IFNγ → ZEB1 → Immune	[58]
PITX2	1.Dev	PITX2	Down	Apo	p38 → PITX2 → Development → Apoptosis	[63, 64]
STAT1	6. Lat	STAT1	Down	-	STAT1 → Development → Immune	[57, 69]
SOX9	7. Unk	SOX9	Down	Apo, Dev	p38 → SOX9 → Apoptosis	[59–62]

TFs changing their concentration levels significantly between wild-type and knock-out mice. “Level”, the changes of TF concentration level in the absence of p38α. Other column headings and abbreviations are the same as those in Table 1

The boldface ones are the main node in Fig. 5

Transcriptional regulatory network for p38α

Gene Ontology (GO) analysis on the target genes of the TFs with strongly changed activity showed enrichment for three GO terms and provided insight into the functional role of the TFs (see Methods, Tables 1 and 2, and Additional file 1: Supplementary Text and Figure S1–S3). The three GO terms are the regulations of the apoptosis (programmed cell death), the downward spiral of the developmental process, and the immune system development. The JNK-c-Jun pathway stimulates the apoptosis, and the I-kB kinase/NF-kB cascade acts as a suppressor of the JNK-c-Jun pathway [41]. Inhibition of p38α MAPK retards another JNK-c-Jun pathway inhibitor NF-kB cascade, but promotes JNK-c-Jun pathway which induces the apoptosis by expressing the Bcl2 protein family [20, 42]. On the other hand, developmental process related genes are down-regulated in the p38α knock-out mice. The study of p38α MAPK [37] reported that the p38α knock-out mice die within days after birth. We do not have enough gene expression profiling data (either other time points in the embryonic period or postnatal period) to investigate TFAs in whole developmental process of the p38α knock-out mice; we cannot confirm but suppose that it might be the reason of the death of the knock-out mice. Further, the genes interact with the TFs which are reported as crucial TFs in the developmental process and the immune system development. Our results are therefore in broad accordance with the experimentally validated results, so it confirmed that our pipeline produces reliable results.

Combining the data in Tables 1 and 2 with those obtained from the literature it is possible to build a putative model for the effects of p38α knock-out (Fig. 5). This figure shows TFs with a strong response in our analysis as nodes, with links that demonstrate regulatory interactions between them. The TF network therefore comprehensively shows the biological consequences of p38α knock-out at transcriptional level.