Identification of interacting transcription factors regulating tissue gene expression in human

Overall analysis procedures

The overall analysis procedures are shown in Additional File 1 and Figure 1. We employed our previously developed function conservation approach [22] to first search for TF pairs using the top 300 expressed genes (referred to tissue-expressed genes) [18] in each of the 79 human tissues from the GNF Atlas2 gene expression database (gnfAtlas2) [23] and their corresponding mouse orthologous genes (Additional File 1). We also utilized promoter sequences from 1018 human housekeeping genes [24] and their mouse orthologous genes to predict interacting TFs playing ubiquitous roles in different tissues (see Methods).

We next filtered out the TF pairs in a particular tissue common to those from housekeeping genes (Figure 1a). The remaining TF pairs (referred to tissue TF pairs) in each tissue were more tissue-specific. The rationale for this filtering is that in each tissue some of the interacting TFs play general roles, since all tissues possess common mechanisms to control the fundamental biological processes. To find interacting TFs controlling tissue gene expression in a broad spatial and temporal manner, we extended the analysis to identify tissue TF pairs common to the same type of tissues (referred to tissue-type TF pairs) as well as interactions of 3 TFs. For the former, we looked for common tissue TF pairs in at least 50% tissues of the same type (Figure 1a). We also built TF-TF interaction networks by joining 2 or more tissue-type TF pairs with one shared TF between TF pairs in the same tissue. TF-TF interaction networks with the same topology in at least 2 tissues from the same tissue type were defined as "tissue-type TF-TF interaction networks" (Figure 1c). Finally, a two-step analysis of TFBS conservation and enrichment of overlapping TF target orthologous genes was performed to predict interactions of 3 TFs (Figure 1b).

Identification of tissue TF pairs

Overall, we identified 2549 tissue and 803 tissue-unique TF pairs for the 79 human tissues. These results indicate that tissue gene expression is regulated by large sets of interacting TFs. Furthermore, the relative small number of tissue-unique TF pairs out of all tissue TF pairs suggests that identical tissue TF pairs in different tissues may play different functional roles, which prompted us to investigate their biological function. For this purpose, we used Gene2go http://www.ncbi.nlm.nih.gov/ to annotate human genes whose promoters contained the target TFBS pairs, as TFs control cellular biological processes via transcriptional regulation of groups of genes with similar functions. Significant (q-value < 0.1) biological processes for tissue TF pairs were obtained by comparing the number of TF target genes involved in a particular biological process to the number of genes for the same biological process in the whole human genome (Fisher's exact test; p = 2.4 × 10^-4 to 8 × 10^-28). All tissue and tissue-unique TF pairs as well as their potential biological functions are listed in Additional File 2.

Evaluation by known interacting TFs

Although the function conservation approach has been proven to be a successful means for predicting interacting TFs [22], we sought to assess the validity of the identified tissue TF pairs. We first used TRANSCompel^® 10.4 [25] to determine if known interacting TFs were statistically enriched in the predicted tissue TF pairs. The TRANSCompel^® database contains 180 experimentally proven composite elements of two or more binding sites which were previously identified by individual wet lab studies from others. Of the180 composite elements, 105 were mapped to the 23,005 (214*215/2) possible combinations of 2 TFs from the 214 non-redundant position weight matrices (PWMs). We first investigated the statistical significance for the occurrence of known interacting TFs in both predicted TF pairs (before filtering) and tissue TF pairs in each of the 79 tissues. Figure 2 shows that known interacting TFs display enrichment more in the tissue TF pairs then in predicted TF pairs for both the number of tissues (37 vs. 9) and the degree of enrichment (Binomial test: p = 3.2 × 10^-2 to 6.8 × 10^-6 vs. p = 4 × 10^-2 to 3.4 × 10^-4). We also computed the occurrence of known interacting TFs in all predicted TFs pairs (before filtering) and all tissue TF pairs from the 79 tissues. We found that 40 (38.1%) of the 105 known interacting TFs were in both predicted TF pairs (Binomial test; p = 6.4 × 10^-9) and tissue TF pairs (Binomial test; p = 5.4 × 10^-11).

To further verify our prediction, we next compared the tissue TF pairs to known tissue-specific TFs from liver, for which the cis-regulatory systems for both individual TF binding and synergistic actions have been thoroughly studied [3, 4, 25, 26]. These studies found 40 liver-specific single TFs and 27 liver-specific interacting TFs. We first computed the tissue TF pairs whose two TFs were all liver-specific based on the 820 (40*41/2) possible combinations of 2 TFs from the 40 liver-specific single TFs. Out of 162 tissue TF pairs from liver tissue, we were able to obtain 30 (Binomial test; p = 2.3 × 10^-14) where both TFs were liver-specific. For the 27 liver-specific known interacting TFs, we found 8 (30%) in both the predicted TF pairs (Binomial test; p = 3.6 × 10^-9) and tissue TF pairs (Binomial test; p = 2.9 × 10^-11) from liver tissue. These include HNF4ALPHA:HNF4ALPHA, NF1:COUP_DR1, CEBPGAMMA:HNF4ALPHA, CEBPA:HNF3B, HNF3:HNF4ALPHA, HNF3:PPARA, CEPBA:GATA4, and HNF1:OCT1. All of them are key elements in liver specific transcriptional regulation. GO enrichment analyses indicated that genes whose promoters contained the predicted liver-specific TFBS pairs were mainly involved in liver specific functions [27, 28], including oxidation reduction, acute-phase response, gluconegnesis, and lipoprotein & lipid metabolic processes (Additional File 2). Further analysis of the binding sites on the promoter sequence of individual genes indicated that we were able to reliably identify interacting TFs similar to those previously reported. One of the examples was the APOA1 gene, which was well-characterized to be synergistically bound by HNF3 and HNF4 [8]. Our prediction was able to successfully identify the HNF3:HNF4 binding sites on its promoter. A closer examination shows that our predicted HNF3 and HNF4 binding sites for the APOA1 gene are exactly those experimentally proven, liver tissue-specific HNF3 and HNF4 binding site combinations [8], which are highly conserved between human and mouse in regards to both nucleotide sequence and spacing between each binding site (Figure 3).

It is important to note that the 79 human tissues represent only part of the temporal and spatial conditions from which the 105 known interacting TFs were discovered, and therefore it is unlikely to have all known interacting TFs in our predicted list. Nevertheless, our results indicate that the use of function conservation approach and tissue-expressed genes was able to reliably identify to a great extent known interacting TFs, thus presenting very strong evidence for the validity of the identified tissue TF pairs. These results also indicate that filtering the TF pairs of housekeeping genes from those in each tissue is an important step to eliminate TFs playing a ubiquitous role, thereby the resulting TF pairs are more tissue-specific.

Identification of tissue-type TF pairs

One of the goals of this study is to find interacting TFs controlling gene expression in a broad spatial and temporal manner such as interacting TFs common to the same type of tissues. This can be achieved by searching tissue TF pairs common across all tissues of the same type such as the 7 muscle tissues. However, the use of all tissues may reduce the power for tissue-type TF pair identification, since the contents of tissue TF pairs and even the function of a common tissue TF pair may be different between tissues of the same type. Therefore, we sought to first classify tissues into smaller but more closely related groups based on tissue TF pairs, from which representative tissues for the same tissue type could be obtained. Accordingly, we used hierarchical clustering to group tissues, as no a priori knowledge was available for the number of groups for each tissue type. The results are shown in Figure 4, where tissues of the same type are generally grouped together such as testis, liver, pancreas, and brain. There are, however, exceptions for other tissue types which are grouped into either distinct groups or into groups with other types of tissues such as muscle and immune systems. While the muscle tissues are classified into two distinct groups, of which one contains skeletal muscle, heart, and tongue and the other contains smooth muscle and cardiac myocytes, tissues for immune systems are classified into a few groups, one of which displays tighten link with a few cancer tissues.

We extended our analysis to investigate conservation for tissue TF pairs between tissues of the two muscle groups. We computed overlap for both tissue TF pairs and their biological functions between tissues using hypergeometric distribution. We found little or no overlap for both tissue TF pairs and their functions among tissues between these two groups, which was especially true for the function of tissue TF pairs (data not shown). On the other hand, both tissue TF pairs and their functions showed significant overlap between tissues within the same group (Figure 5). These results not only demonstrate the validity of our tissue classification but also indicated that tissues from the same type (here 5 muscle tissues) may have great difference in both the contents of tissue TF pairs and TF functional roles.

Reconstruction of tissue-type TF-TF interaction networks

In an effort to reveal TF relationships in controlling tissue gene expression, we performed analysis to reconstruct TF-TF interaction networks. Using tissue-type TF pairs, we first looked for those with one shared TF between each other in the same tissue, from which TF-TF interaction networks were built by joining 2 or more TF pairs (Figure 1c). TF-TF interaction networks with the same topology in at least 2 tissues from the same tissue-type group were then selected as tissue-type TF-TF interaction networks, which are multi-input network motifs consisting of at least 3 TFs that bind to a set of gene promoters. A total of 84 tissue-type TF-TF interaction networks were identified for the 11 tissue-type groups, ranging from 1 for immune/cancer to 22 for testis (Additional File 4). Sixty two of these tissue-type TF-TF interaction networks have a linear relationship between TFs with 1 to 4 internal TFs (i.e. TF connecting to 2 other TFs), indicating that the majority of the TF-TF regulatory networks have simple TF relationships for controlling tissue gene expression. Figure 6a shows a multi-input network motif from liver tissues, in which FOXJ2, HNF1, and TTF1 regulate 6 genes in a combinatorial manner by either 2 or 3 TFs. The remaining 22 tissue-type TF-TF interaction networks display more complex interacting structures with some of the internal TFs connecting to 3 or more TFs.

Unlike the tissue TF pairs, we did not find any common tissue-type TF-TF interaction networks between different tissue types. In light of this, we performed a search to see if any single TFs played central roles in controlling tissue gene expression across different tissues, and looked for internal TFs in multiple tissue-type TF-TF interaction networks. To this end, we found that AP2, PPARA, PAX4, FAC1, ZIC3, and SPZ1 served as internal TFs in 8, 8, 8, 6, 5, and 4 tissue-type TF-TF interaction networks, respectively, suggesting their role as central hubs in tissue-type TF-TF interaction networks. Whereas FAC1 acts as the internal TF in 6 tissue-type TF-TF interaction networks from immune systems and cancer, SPZ1 mainly serves as the internal TF in tissue-type TF-TF interaction networks from testis, and the rest in 5 to 6 tissue-type TF-TF interaction networks from different tissue types. These results indicate that FAC1, when serving as the internal TF, is restricted to the two related tissue types, and that SPZ1, a bHLH-Zip protein, has an important role in testis [29, 30]. The rest have more diversified roles for coordinating network TFs in controlling tissue gene expression.

It is interesting to note that no single TFs serve as the central hub for tissue-type TF-TF interaction networks from liver tissue. However, we observed that 6 of 7 tissue-type TF-TF interaction networks had at least one known liver-specific TF serving as the internal TF as shown in Figure 6b. To investigate if this distribution pattern of liver-specific TFs in the TF-TF interaction networks had any biological meaning, we randomly sampled TFs from the 214 PWMs to build TF-TF networks, each having the same size and order as the real TF-TF interaction networks. The simulated TF-TF networks were then compared to tissue-type TF-TF interaction networks to estimate the statistical significance for the distribution of liver-specific TFs. The results indicate that known liver-specific TFs were significantly enriched as internal TFs for these 7 tissue-type TF-TF interaction networks (bootstrap analysis; p < 10^-20). By contrast, the total number of liver-specific TFs in these 7 tissue-type TF-TF interaction networks was not enriched (bootstrap analysis; p = 0.11). These results suggest that liver-specific TFs, other than initiating liver-specific transcriptional event, may play an important role in recruiting non-liver-specific TFs to the tissue-type TF-TF interaction network, thus offering the potential for coordinating and controlling gene expression across a variety of conditions.

Prediction of multiple interacting TFs

As a first step to elucidate cis-regulatory modules involved in tissue gene regulation, we extended our analysis to the interactions of 3 TFs (named as multiple interacting TFs). Using tissue TF pairs from each of the 79 tissues, we performed a two-step analysis of TFBS conservation and enrichment of overlapping orthologous genes between human and mouse (see Methods). Although it is likely that multiple interacting TFs may be under estimation by the use of tissue TF pairs instead of all predicted TF pairs, the predicted multiple interacting TFs are tissue-specific. Therefore, these predictions most likely represent cis-regulatory modules involved in tissue gene regulation. To this end, we identified 1735 unique interactions of 3 TFs for the 79 human tissues, ranging from 9 multiple interacting TFs for testis interstitial to 72 multiple interacting TFs for caudate nucleus (Additional File 5).

The validity of these predicted multiple interacting TFs was assessed by using liver-specific single TFs from TRANSCFAC11.4 [25], as few known cis-regulatory modules were available. We performed analysis to see if known liver-specific TFs were statistically enriched in 30 predicted multiple interacting TFs from liver tissue. We found 4 of them (bootstrap analysis; p < 10^-3) whose 3 TFs were all liver-specific, 18 (bootstrap analysis; p < 10^-8) with at least 2 liver-specific TFs, and 28 (bootstrap analysis; p < 10^-5) with at least 1 liver-specific TF. These results provide evidence for the enrichment of liver-specific TFs in the predicted multiple interacting TFs, which in turn demonstrated the validity of the prediction.

We next searched for all predicted multiple interacting TFs and their potential functions that are common between tissues. The results indicated that, although common multiple interacting TFs existed between most tissues, the highest overlap was within brain tissues and between brain and gland tissues. By contrast, there was little overlap for the functions of multiple interacting TFs, except within brain and cancer and between these 2 tissue types (Additional File 6). The latter is especially interesting to us, as cancer cells have a global effect on immune systems, which in turn control and shape developing cancer [31]. Six multiple interacting TFs were found to have common functions between immune systems and cancer tissues, including CEBPGAMMA:NKX25:PLZF, CEBPGAMMA:PAX4:PLZF, CP2:NFY:PAX4, FOXJ2:PAX4:POU3F2, CEBPGAMMA:PAX4:PLZF, and FOXJ2:HNF3:PAX4. These results revealed not only the common mechanisms for transcriptional regulation but also the common functional role of multiple interacting TFs between cancer and immune systems, including cell cycle, cell division, DNA replication, mitosis, phosphoinositide-mediated signaling, and immune response. These findings therefore provide new insight into the molecular interplay between cancer and immune systems.

Promoter sequences for housekeeping and tissue-expressed genes

The GNF Atlas2 gene expression database (gnfAtlas2) [34], which contains gene expression data from 79 human tissues, was used for the selection of genes. Based on the report of Pennacchio et al[18] we selected in each tissue the top 300 expressed genes (referred to tissue-expressed genes), which have been used and proven to successfully predict tissue-specific enhancers. Housekeeping genes are the 1018 genes defined by Farre et al[24]. Redundant genes in each group were first removed. Although regulatory elements can exist anywhere in the genome, they are more concentrated around the transcriptional start sites [35]. To reduce false predictions we focused on the proximal promoters which have been proven to successfully predict tissue-specific regulatory elements [21, 36]. It is however worthy to note that the use of 1 kb promoter sequences has limitation for the prediction of tissue TF pairs when compared to the experimental approaches such as ChIP-chip experiment, in which TF pairs can be detected anywhere in the genome. Considering no benchmark promoter sequence dataset is currently available for computational prediction of functional TF pairs, the use of 1 kb promoter sequences and our computational approach nevertheless provide an alternative way for tissue interacting TF discovery. Promoter sequences within 1 kb upstream of transcriptional starting sites for both human and corresponding mouse orthologous genes were extracted from the UCSC genome browser (hg18 March 2006 assembly, mm9 July 2007 assembly). Orthologous genes with promoter sequences from both human and mouse were selected for further analysis. This procedure resulted in 208 to 278 orthologous promoter sequences for tissue-expressed genes and 986 orthologous promoter sequences for housekeeping genes.

Prediction of TF pairs and tissue TF pairs

The procedures for predicting TF pair are basically the same as previously described [22] (Additional File 1). Briefly, background sequences were created by shuffling the DNA sequences within each promoter by either mixing completely or keeping dinucleotides together. These background sequences are preferable to using intergenic sequences which usually are AT-rich or exonic sequences whose nucleotide distributions tend to be biased, when compared to the test promoter sequences. The resulting shuffled sequences from human and mouse, together with the original promoter sequences, were employed for TFBS detection using the Match^® program [37], for which the profile parameter was set to "minimize the sum of false positives and negatives", and the 214 non-redundant vertebrate PWMs from the professional TRANSFAC11.4 database [25]. To detect enriched TF pairs out of 23,005 (214*215/2) possible combinations of 2 TFs, distance constraints were first applied for the selection of co-occurring TFBSs with a defined maximum distance between 2 TFBSs. A total of 10 distances were defined, ranging from the smallest 20 bp to the largest 200 bp with a 20 bp increment. The assumption behind the distance constraint is that functional TFBS pairs are more distance-restricted than random co-occurrence of TFBSs [19, 38]. This is true not only in human, for which we found that functional TF pairs were enriched within 200 bp distance ranges [22], but also in Drosophila, in which short-range linkages (< 50 bp) between TFs was overrepresented but mid-range distances (100-500 bp) between TFs was depleted [39]. Enrichment of TFBS pairs for each distance constraint was achieved by computing the ratio of counts for a particular TFBS pair in real promoter sequences vs. the counts of the same TFBS pair in background sequences. To reduce noise while keeping as many as TFBS pairs for the integration of function conservation analysis described below, TFBS pairs with ratio > 1 in more than 5 distance constraints were selected.

A two-step analysis procedure was employed to compute the enrichment of overlapping orthologous genes for a particular TFBS pair. First, a cutoff threshold of at least 10% overlapping orthologous genes between mouse and human was set up for selecting genes whose promoters contained the TFBS pair. The enrichment of overlapping orthologous genes was then estimated by computing the ratio of overlapping orthologous genes from real promoter sequences against those from shuffled sequences. This analysis was performed for each distance constraint. The integration of function conservation for each TF pair was achieved by estimating the correlation (Pearson correlation coefficients) between the 10 enriched TFBS pairs and 10 corresponding enriched overlapping orthologous genes from the same distance constraint. Permutation tests were employed to estimate the statistical significance of correlation by randomly matching the 10 TFBS pair ratios with the 10 overlapping orthologous gene ratios. For multiple test correction, a cutoff threshold of q-value < 0.05 was applied. TF pairs are those passing the cutoff and common between human and mouse.

We next filtered TF pairs of housekeeping genes from those in each tissue (Figure 1a). This was done by removing TF pairs in a particular tissue common to those from housekeeping genes. The remaining TF pairs in each tissue were more tissue-specific, and therefore, were defined as tissue TF pairs. Similar results were obtained from using background sequences of either completely mixed nucleotides or keeping dinucleotides together or completely mixing nucleotides. The results from completely mixed nucleotides were used.

Clustering analysis

To group tissues based on their tissue TF pairs, a 2549 (tissue TF pairs) × 79 (tissues) matrix with binary numbers was first built for all tissue TF pairs in the 79 human tissues. The presence of a tissue TF pair in the matrix was labeled with 1 and the absence was labeled with 0. A distance matrix was then built using the "binary" method, and hierarchical clustering was subsequently performed using the "complete" agglomeration method. All analysis was performed using the R statistical package [40].

Predicting multiple interacting TFs

A two-step analysis of TFBS conservation and enrichment of overlapping orthologous genes was performed to predict interactions of 3 TFs. For TFBS conservation, the identified tissue TFBS pairs were first used to construct all possible 3-TFBS combinations by searching paired tissue TFBS pairs with one shared TFBS between each other on exactly the same location of a gene's promoter (Figure 1b) in a particular tissue. Orthologous gene pairs containing conserved 3-TFBS combination between human and mouse were then selected. Conserved 3-TFBS combinations are those whose 3 TFBSs have the same order and orientation on the promoter sequences between human and mouse orthologous genes. For enrichment of overlapping orthologous genes in a tissue, however, multiple interacting TFs from different orthologous gene pairs were considered to be the same as long as they contained the same 3 TFs. Enriched multiple interacting TFs are those with 3-TFBS combinations occurring on at least 10 orthologous gene promoters and with their target orthologous genes displaying significant overlap between human and mouse (p = 3 × 10^-2 to < 10^-36 and q < 0.05).

Statistical methods for enrichment analyses

Two main statistical methods were employed for estimating the significance of enrichment in this study. For validating predicted TF pairs and tissue TF pairs by known interacting TFs, the binomial distribution probability, as shown below, was used to determine if known interacting TFs were present more often in the predicted TF pairs or tissue TF pairs than in a randomly selected group from a given list of TFs.

For example, in the case of estimating the statistical significance of known liver-specific interacting TFs in our predicted tissue TF pairs from liver tissue, the n is the number of known liver-specific interacting TFs in the predicted tissue TF pairs from this study; N the number of tissue TF pairs from liver tissue; and p _f the background probability of liver-specific TF pairs in all possible combinations of 2 TFs from 214 PWMs.

The statistical significance was computed using the hypergeometric distribution to estimate (1) the enrichment of overlapping tissue TF pairs and their overlapping functions between muscle tissues, and (2) overlapping orthologous genes in predicting multiple interacting TFs.

In the case of overlapping orthologous genes in predicting multiple interacting TFs, for example, c is the number of orthologous gene pairs containing conserved 3-TFBS combination between human and mouse; N the number of tissue-expressed genes for a particular tissue; S₁ and S₂ are the numbers of tissue-expressed genes with 3-TFBS combinations corresponding to those in c for human and mouse, respectively. The resulting p-value is the probability of observing c or more orthologous gene pairs containing conserved 3-TFBS combination from two sets of size S₁ and S₂ drawn from a set of N tissue-expressed genes.