Defining transcriptional networks through integrative modeling of mRNA expression and transcription factor binding data

Abstract

Background

Functional genomics studies are yielding information about regulatory processes in the cell at an unprecedented scale. In the yeast S. cerevisiae, DNA microarrays have not only been used to measure the mRNA abundance for all genes under a variety of conditions but also to determine the occupancy of all promoter regions by a large number of transcription factors. The challenge is to extract useful information about the global regulatory network from these data.

Results

We present MA-Networker, an algorithm that combines microarray data for mRNA expression and transcription factor occupancy to define the regulatory network of the cell. Multivariate regression analysis is used to infer the activity of each transcription factor, and the correlation across different conditions between this activity and the mRNA expression of a gene is interpreted as regulatory coupling strength. Applying our method to S. cerevisiae, we find that, on average, 58% of the genes whose promoter region is bound by a transcription factor are true regulatory targets. These results are validated by an analysis of enrichment for functional annotation, response for transcription factor deletion, and over-representation of cis-regulatory motifs. We are able to assign directionality to transcription factors that control divergently transcribed genes sharing the same promoter region. Finally, we identify an intrinsic limitation of transcription factor deletion experiments related to the combinatorial nature of transcriptional control, to which our approach provides an alternative.

Conclusion

Our reliable classification of ChIP positives into functional and non-functional TF targets based on their expression pattern across a wide range of conditions provides a starting point for identifying the unknown sequence features in non-coding DNA that directly or indirectly determine the context dependence of transcription factor action. Complete analysis results are available for browsing or download at http://bussemaker.bio.columbia.edu/papers/MA-Networker/.

Background

For various organisms, DNA microarrays have been used to measure the mRNA abundance for essentially all protein-coding genes in the genome under a large number of conditions [1, 2]. Microarray technology can also be combined with chromatin-immunoprecipitation (ChIP) or chromatin profiling (DamID) to quantify the occupancy of upstream non-coding regions by transcription factors or other chromatin-associated proteins [3–9]. In the budding yeast Saccharomyces cerevisiae, ChIP has been used to globally map the binding sites of over a hundred transcription factors [4]. Moreover, mRNA expression data for over a thousand conditions has been published. The challenge is to find new ways to extract knowledge about the regulatory mechanisms that govern the cell by combining these different types of data [10–14].

Initiation of transcription in eukaryotes is a complicated process that depends on the binding of transcription factors (TFs) and chromatin-modifying enzymes to the promoter region as well as the recruitment of the RNA Polymerase II complex to the transcription start site. Transcriptional control is combinatorial, and cooperative binding of multiple factors on the same promoter region and/or cooperative recruitment of the Pol II complex is often required for transcriptional activation [15]. Occupancy of the promoter region of a gene by a transcription factor is thus a necessary but not a sufficient condition for the gene to be controlled by it. As a consequence, genome-wide transcription factor binding patterns measured using ChIP or DamID microarray experiments alone can only indicate the potential for a gene to be regulated by a given TF. Independent information will be required to establish that the gene is indeed a functional target of the factor.

Deletion or over-expression of a transcription factor, combined with genomewide microarray profiling of the difference in expression between mutant and wild type, is also widely used to infer regulatory interactions. However, drastic perturbation of the genetic network outside the physiologically relevant range may lead to false target prediction, or the mutant strain may simply not be viable. Moreover, direct and indirect targets of the factor cannot be distinguished using this approach.

When mRNA abundances for all genes are compared between two experimental conditions in a microarray experiment, the observed differential expression pattern is usually a superposition of responses of various pathways, mediated by signaling cascades that end at the level of transcription factors. It has previously been shown that these changes in TF activity can be quantitatively inferred by performing multivariate regression analysis on the expression log-ratios from a single microarray experiment [16–19]. Transcription factors are implicitly represented by a consensus motif for their DNA binding sites, and the regression coefficients estimate the changes in TF activity.

In the present study we will instead use ChIP occupancy log-ratios as predictors for expression. No sequence information will be used, and it is therefore not necessary that a DNA consensus motif be known for the TFs. Again, multivariate regression analysis of a single genomewide set of mRNA log-ratios on the genomewide binding profiles of a large number of TFs for which ChIP data is available can be used to quantify to what extent each transcription factor is responsible for the observed changes in mRNA expression.

When the regression procedure is performed in parallel on a large library of expression data, the inferred changes in TF activity for each comparison can be combined into a transcription factor activity profile (TFAP) for each transcription factor ("Step 1" in Fig. 1a). Each TFAP represents a highly specific regulatory signature, as is shown for three transcription factors in Fig. 1b: The activity of the G2 phase related factor Ndd1p oscillates during the cell cycle (blue), but shows little or no response to nutrient depletion and other stress conditions (red), or changes in alpha pheromone concentration (green). Complementary behavior is seen for the TCA cycle regulator Hap4p and the mating-related factor Ste12p.

Figure 1

(a) Overview of our method for determining regulatory coupling strengths between transcription factors and their putative target genes. Inputs are (i) a library of microarray expression data for a large number of conditions and (ii) genomewide (ChIP) occupancy data for one or more transcription factors. In the first step of our algorithm, a matrix of transcription factor activities is inferred by using regression analysis to explain the mRNA expression pattern under each condition in terms of the ChIP data for each transcription factor. In the second step, a matrix of regulatory coupling strengths is determined by computing the correlation between each transcription factor activity profile (TFAP) and the mRNA expression profile of each gene. (b) Examples of transcription factor activity profiles. The activity profiles of three transcription factors (Hap4, Ndd1, Ste12) are shown across stress response, pheromone response, and cell cycle [28–30]. Significant changes in activity of the TCA cycle regulator Hap4p occur mostly in metabolic stress conditions, while changes in the activity of the cell cycle regulator Ndd1p and the pheromone-dependent regulator Ste12p are associated with the cell cycle and signal transduction experiments, respectively. (c) Examples of scatter plots of ChIP binding log-ratio versus coupling factor. In the scatter plots, black dots denote unbound (B-) genes, red dots denote bound and coupled genes (B+/C+), while green dots denote genes that are bound but not coupled (B+/C-). A threshold of P = 10^-3 was used to determine significant binding as described in Lee et al. [4]. A threshold for coupling was determined by requiring a false discovery rate of 5%, as described in Methods.

One expects the mRNA expression profile of a gene regulated by a specific transcription factor to be similar to the TFAP of that factor. We therefore investigated whether the linear correlation across the experiment library between a TFAP and the mRNA expression profile of a gene whose promoter is bound by the factor could be interpreted as a regulatory coupling strength and used to improve the specificity of target prediction. To this end, we constructed a matrix of regulatory coupling strengths between all transcription factors and all genes ("Step 2" in Fig. 1a). When this information is combined with the original ChIP data for a given TF, the ChIP log-ratio and coupling strength for each gene can be shown simultaneously in a 2D scatter plot (Fig. 1c). The fact that each gene has two parameters associated with it allows a more sophisticated classification than is possible based on ChIP alone. We first defined a set "B+" of genes that are significantly bound by a TF (we required the P-value reported by Lee et al. to be smaller than 10^-3) [4]. Next, we then partitioned the "B+" gene set into two subsets "B+/C+" and "B+/C-" based on whether or not their mRNA expression profile was significantly correlated with the TFAP (Pearson correlation, 5% false discovery rate). Our hypothesis is that the B+/C+ genes (shown in red in Fig. 1c) are the functional direct targets of the factor, while the binding to the promoter region of the B+/C- genes (shown in green) is non-functional.

Results and discussion

Only a subset of genes bound by each TF is controlled by it

We have focused on the yeast Saccharomyces cerevisiae, for which a wealth of functional genomics data is available. We compiled a library of ~750 expression patterns from various sources, and combined it with the genome-wide promoter occupancies in mid-log phase and rich medium for 113 TFs as measured by Lee et al. [4]. It was determined that 37 transcription factor occupancy patterns out of the full set of 113 are significant predictors of mRNA expression for one or more experiments in our library (see Methods). Note that the library of expression data we used obviously does not cover all possible experimental conditions and our method is therefore likely to underestimate the number of transcription factors that are present in the nucleus under the conditions used by Lee et al. [4]. For each of the 37 factors selected for further analysis, a transcription factor activity profile (TFAP) was computed ("Step 1" in Fig. 1a) and the TF-target coupling strength was determined ("Step 2" in Fig. 1a). The number of genes in each group (B-, B+/C+, and B+/C-) is listed in Table 1 for the 37 transcription factors analyzed. On average 58% of significantly bound genes are classified as significantly coupling genes. Activity profiles and B+/C+ target predictions for all 37 factors are available on the website supporting this paper.

Table 1 Classification of genes according to ChIP data and inferred regulatory coupling.

Enrichment for specific functional categories

Several analyses were performed to validate our results. First, we established that B+/C+ genes are significantly enriched for specific Gene Ontology (GO) categories (hypergeometric distribution; 5% false discovery rate) [20]. This result is not surprising, as we would already expect the set B+ of ChIP positives per se to be enriched for roughly the same functional categories. By contrast, for almost all TFs analyzed we found no significant enrichment for any GO category in the set of non-coupling (B+/C-) genes (Fig. 2, see supplementary website for details). This result is very significant because it suggests that our criterion for distinguishing functional from non-functional TF targets based on regulatory coupling is accurate: There seems to have been no evolutionary pressure on the set of B+/C- genes.

Figure 2

Enrichment for functional annotation. The number of significantly over-represented Gene Ontology (GO) categories in the group B+/C+ of genes that couple to transcription factor activity (red) and the group B+/C- of genes that do not couple (green) for each of the 37 transcription factors analyzed. No significant enrichment for any GO category was found in most B+/C- gene groups, supporting the hypothesis that only the coupling B+/C+ genes are functional targets.

Transcriptional response to transcription factor deletion

Next, we analyzed the expression response to transcription factor deletion for the B+/C+ and B+/C- genes. The mean and the standard deviation of the gene expression log-ratio between mutant and wild type as obtained in Hughes et al. were calculated for all genes in the genome, as well as for the B+/C+ and B+/C- groups [21]. A sample t-test was performed to determine the significance of the change in expression. The B+/C+ genes show a significant change in mRNA expression for the 7 transcription factors for which deletion and ChIP data are both available. By contrast, the response of the B+/C- genes is insignificant for most transcription factors (Table 2).

Table 2 Analysis of transcriptional response to transcription factor deletion.

Enrichment of promoter regions for consensus binding motifs

In a third analysis to validate our results, we tested for over-representation of 4 different cell cycle related DNA consensus motifs (MCB, SCB, Swi5p and SFF) in the upstream regions of 7 cell cycle related TFs (Ace2, Fkh2, Mbp1, Mcm1, Ndd1, Swi4 and Swi5). The binomial distribution was used to score motif enrichment in the B+/C+ and B+/C- gene sets for each of the transcription factors. We found certain DNA motifs to be significantly over-represented in B+/C+ genes for one or more TFs, but dramatically less so in B+/C- genes (Table 3). Finally, we defined B+/C+ groups using duplicate ChIP experiments for 7 cell cycle regulators performed by Simon et al. and found the overlap with the B+/C+ genes for the data of Lee et al. to be 85% on average [4, 22].

Table 3 Over-representation of four cell cycle related motifs.

Assigning directionality to divergently transcribed promoters

Taken together, the results mentioned above convincingly demonstrate that the use of a coupling factor threshold as a novel additional criterion leads to significantly improved specificity in the prediction of functional TF targets. The biological implications of our analysis are highlighted in the case of divergently transcribed genes that share a common promoter region, represented as a single microarray probe. There are 1592 such probes out of the total 4532 probes in the ChIP experiments of Lee et al. [4]. When the ChIP data indicate that a TF binds to the intergenic region, nothing can be said about whether it regulates one of the genes or both based on that information alone. By contrast, our regulatory coupling analysis naturally allows us to distinguish between these different scenarios and make precise statements about which genes are controlled by each of the factors that occupy the promoter region (see Fig. 3). Both uni- and bi-directional control by TFs is observed. Indeed, we found the functional annotation of the protein encoded by the coupled targets to be consistent with what was known about the function of the bound TF in most cases analyzed [20].

Figure 3

Assigning directionality to divergently transcribed promoters. For pairs of divergently transcribed genes sharing a single promoter region occupied by one or more transcription factors, our method can be used to determine which gene is regulated by which factor. In the diagrams, genes are represented as squares with arrows showing the transcription direction; transcription factors are shown as ovals. The numbers shown are significance scores for the coupling between the transcription factor and the gene, equal to the negative 10-based logarithm of the P-value. Significant regulatory relationships are shown as arrows to colored boxes. In (A), the cell cycle transcription factor Mbp1p regulates the recombinase RAD51 but not the endopeptidase PUP3, while in (B), the putative rRNA processing regulator Fhl1p regulates the ribosomal subunit RPL40B but not the protein kinase MLP1. In the scenario illustrated in (C), both Nrg1p and Hap6p bind to the intergenic upstream region of FSP2 and HXT9. The coupling analysis shows that Nrg1p in this case works bi-directionally and regulates both genes, while Hap6p regulates neither gene.

Revealing intrinsic limitations of TF deletion experiments

Since TF occupancy data from ChIP experiments can be used to separate direct from indirect targets among the genes that respond to TF deletion, combining ChIP data with deletion data can potentially achieve the same goal as our more sophisticated analysis. Keeping Occam's Razor in mind, it is therefore important to investigate to what extent mRNA expression log-ratios from a TF deletion mutant vs. wild type experiment can replace our regulatory coupling strength on the vertical axis in Fig. 1b. We defined sets B+/D+ and B+/D- for each of the seven TFs for which data are available in Hughes et al., the D+ genes being those that show a response to the TF deletion at P < 10^-2, using the P-value provided by those authors [21]. In the regulatory coupling analysis described above, we found no significant enrichment for specific GO categories in the group B+/C- of genes that are bound but not coupled. In the present case however, similar levels of enrichment for GO categories are found for B+/D+ and B+/D- genes (Table 4). This result indicates that a substantial fraction of the functional direct targets of a typical transcription factor is being missed even if one combines transcription factor deletion data with ChIP data. A plausible explanation for this lack of sensitivity is that each TF deletion experiment, by definition, is performed under a single condition in which not all possible co-factors of the deleted TF may be present. By using a large and heterogeneous library of experimental conditions as input, our method samples most or all co-factor combinations that occur as the context for control by each TF, naturally taking into account the combinatorial nature of transcriptional control.

Table 4 Replacing regulatory coupling strength by response to transcription factor deletion.

Conclusions

Our results underscore the unique added value of ChIP data such as that of Lee et al. when it is used in combination with a library of mRNA expression data [4]. We found that roughly half of the transcription factor targets predicted by ChIP are nonfunctional. Although some of these will be false positives of the ChIP technology, especially for TFs that are not present in active form in the nucleus under the conditions used by Lee et al., we believe that our results instead point to interesting biology: TF binding can fail to confer transcription of a nearby gene for a variety of reasons, including competition with nearby activators or repressors, local or global chromatin conformation, or lack of partners for cooperative recruitment of the Pol II complex.

Several works have relied on representing a TF by its mRNA expression profile in order to discover connections between transcription factors and their targets [23–25]. By contrast, our method infers changes in TF activity by analyzing the mRNA levels of putative TF targets. This allows us to analyze regulatory relationships even if the TF is modulated in a purely post-translational manner, e.g. by phosphorylation. The reliable classification of ChIP positives into functional and non-functional TF targets, as it has been presented here, provides a starting point for future research aimed at identifying the unknown sequence features in non-coding DNA that directly or indirectly determine the context dependence of TF action.

Methods

Microarray expression and binding data

A library of 751 genomewide mRNA expression patterns (transcriptomes) was compiled from a variety of sources (see supplementary data for complete references). ChIP data for 113 transcription factors was downloaded from the website accompanying Lee et al. [4]. We used the P-values provided by these authors to determine which genes were significantly bound by each given factor at a confidence level of P < 10^-3. All microarray data used in our analysis was represented as log-ratio base two.

Transcription factor activity profiles

For each separate transcriptome t, we used the following multivariate regression model to infer transcription factor activities for each microarray experiment:

Here E _gt represents the mRNA expression log-ratio of gene g in experiment t, while B _fg represents the ChIP log-ratio for transcription factor f and the promoter region of gene g. The intercept F_0trepresents a baseline expression level, while the regression coefficients F _ft can be interpreted as inferred transcription factor activities. Starting with the full set of 113 transcription factors, we used backward selection to eliminate uninformative transcription factors from our model: First, for each microarray experiment a P-value corresponding to each regression coefficient was determined, based on an F-test [26]. The transcription factors were then sorted based on the smallest P-value among all 751 experiments. In an iterative procedure, the transcription factor with the most insignificant P-value was removed until all factors were significant at a P-value of 0.005/751. Since this analysis in itself is novel and useful, we have made an online ChIP regression tool available at http://bussemaker.bio.columbia.edu/tools/.

Gene-TF coupling factor

For each pair-wise combination of a gene g (represented by its mRNA expression profile E _gt ) and a transcription factor f (represented by the inferred activity profile F _ft ), a regulatory coupling factor was calculated, equal to the Pearson correlation between E _gt and F _ft :

For each value of r, an associated P-value was computed by performing a t-test on t = r[(G-2)/(1-r²)]^1/2. To account for the parallel testing of many TF-target pairs, but at the same time avoid the overly conservative Bonferroni correction, we set a threshold for t by requiring a false discovery rate of 5% [27]. The end result of this procedure is a list of genes that are significantly coupled to a transcription factor. Strictly speaking, to avoid circularity, the coupling of each gene should be evaluated based on a TFAP derived from expression data for all but that gene. However, as the TFAP is derived from the expression profile of all genes bound by the TF, the effect of leaving out one gene is relatively insignificant in practice. Moreover, repeating this procedure for every gene in the genome would be computationally unfeasible.

Enrichment for gene ontology categories

Based on the regulatory coupling analysis described above, the genes bound a given transcription factor (B+) were sorted in two classes, B+/C+ (bound and coupled) and B+/C- (bound but not coupled). These two sets were used as input for further analysis. The cumulative hypergeometric distribution was used to determine whether a set of genes is enriched for one or more Gene Ontology categories [20]. The Bonferroni correction was applied to all P-values to deal with the parallel testing of GO categories. The organism-independent ontology and the gene-association table (version May 2003) for S. cerevisiae were downloaded from http://www.geneontology.org.

Response to transcription factor deletion

Expression data for mutant vs. wild-type comparison for the transcription factors Dig1, Gcn4, Hir2, Mbp1, Swi4, Swi5, and Yap1 were obtained from Hughes et al. [21]. To test whether a given subset of genes responded to TF deletion, a sample t-test was performed, comparing the average expression log-ratio in the subset with the genome-wide distribution of expression changes. To guarantee that this analysis was fair, the respective TF deletion experiments were excluded from the library used to calculate the coupling factors that define the C+ and C- groups.

Enrichment for cell cycle DNA motifs

Four different DNA motifs found as top-scoring motifs by REDUCE and also reported in Spellman et al. were tested for over-representation in the set of B+/C+ and B+/C- genes, respectively, for the 7 cell-cycle related transcription factors within the set of 37 factors analyzed by us [16, 28]. These motifs are: ACGCGT (MCB), CGCGAAA (SCB), AACCAGC (Swi5p) and GTAAACA (SFF). Motifs were counted in non-coding regions up to 600 bp upstream from the ORF start position, and expected counts were based on upstream regions of all genes. No overlapping matches were counted. The cumulative binomial distribution was used to assign a P-value to the enrichment for these motifs.