Abstract

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. ​Fig.1a).1a). Each TFAP represents a highly specific regulatory signature, as is shown for three transcription factors in Fig. ​Fig.1b: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.

Open in a separate window

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. ​Fig.1a).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. ​(Fig.1c).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. ​Fig.1c)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

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. ​(Fig.2,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.

Open in a separate window

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.

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. ​Fig.3).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].

Open in a separate window

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.

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

List of abbreviations

TF: transcription factor

ChIP: chromatin immunoprecipitation

DamID: DNA adenine methyltransferase identification

TFAP: transcription factor activity profile

GO: gene ontology.

Authors' contributions

FG and HJB both contributed to the development of the algorithm and the analysis and presentation of the results. BCF compiled the library of expression data used and contributed tools for functional annotation enrichment analysis. All authors read and approved the final manuscript.

Acknowledgements

References