Two different classes of co-occurring motif pairs found by a novel visualization method in human promoter regions

The transcription of genes is considered to be primarily regulated by transcription factors (TFs). TFs are recruited to interact with specific DNA sequences or motifs in promoters. It is an important and challenging issue to characterize these promoter motifs. It was shown that some motifs (e.g., TATA, SP1, and CREB) are most commonly found within -400 to +100 bp of the transcription start site (TSS) [1]. In other words, positional distribution plots of single motifs found in many promoters show a single peak in frequency near the TSS [1, 2]. In the promoters of higher eukaryotes, such as humans, multiple TFs act in a coordinated way to enable complex patterns of gene expression; thus, it is highly likely that even those binding sites of TFs that do not show clear positional preferences are placed under some constraints, showing a bias such as a frequency peak at a certain distance from the TSS when the co-occurrences of a motif pair are considered. To date, no study has examined the positional relationships among motif pairs and TSS. Furthermore, the search for motifs is usually limited to the 1 kb of sequence immediately upstream of a TSS; however, more than 1 kb of upstream sequence is commonly conserved among vertebrates. A previous study has reported a diverse range of different TSSs for the same genes [3], and distal-promoters are also important for gene expression [4]. Therefore, > 1 kb of sequence should be examined for motif pairs.

In higher eukaryotes, TF binding sites are often organized in clusters called cis-regulatory modules (CRM) [5]. Most of the computational CRM prediction algorithms developed to date rely mainly on the distances between motifs and the phylogenetic conservation of these motifs across different species [6–8]. For example, Blanchette et al. recently performed a genome-wide prediction of CRMs in the human genome [8]; however, further analyses are required to clarify the nature of different motifs, such as which motif pairs are preferentially used in promoters. Yu et al. [9] addressed this issue by identifying motif pairs with biased distance distributions in various sets of human genes that exhibited tissue-specific expression, and Long et al. [10] identified motif pairs over-represented in promoters of immune-response genes. It also appears to be important to gain an understanding of the function of upstream sequences that are not specific to differentially expressed genes. Such overall analyses may reveal the specific features of (a partial set of) house-keeping genes.

We are currently engaged in the construction of two databases (DBTSS [11] and H-InvDB [12]) that contain information on TSSs determined by the mapping of cDNAs with intact 5'-ends onto genome sequences. H-InvDB covers a wide range of genes, but the TSS information is expected to be more accurate in DBTSS. This wealth of data should be useful in obtaining further insights into motif positioning relative to TSSs.

In this study, we developed a method of identifying positional relationships among motif pairs and TSSs in human promoters. We searched from -2000 to +1000 bp relative to each TSS, enabling an examination of possible alternative promoters [13] or distant motifs in enhancers and silencers. We first visualized the positional preferences among pairs of non-overlapping motifs and the representative TSS for each gene. This visualization assisted in the identification of two distinct classes of motif pairs in terms of the positional relationships between the pairs and TSSs. Differences between these two classes were examined with regard to the positional distributions of their constituent motifs, the expression patterns of the genes that harbor the motif pairs, and the protein domains relevant to the transcription factors that bind to the promoter motifs.

We identified TRANSFAC motif pairs that significantly co-occur within human promoter regions. This procedure is dependent on the accuracy of the TSS information. For TSS annotation, we used two databases: H-InvDB and DBTSS. H-InvDB contains more sequence data of transcripts that cover not only protein-coding mRNAs but also Pol II-transcribed non-coding RNAs. Because the promoter sequences are more abundant in H-InvDB, fewer false negatives are expected; it is therefore natural that a greater number of significant motif pairs were found in H-InvDB than in DBTSS. On the other hand, the advantage of using DBTSS lies in the accuracy of its TSS information. Given that we confirmed that the co-occurrence patterns for motif pairs derived from both annotations were similar, we used H-InvDB for down-stream analyses.

The detection ratio of significantly detected motif pairs was 69% (86/127), while we detected 32% of the total possible combination (i.e. 8,454 motif pairs out of 228 × 227/2 motif pairs), which is a promising result considering the previous work by Yu et al. [9]. Its detection ratio (or sensitivity) was 40%, detecting 19% (9,060 motif pairs) of 46,971 combination using tissue specific gene groups. Our detection ratio is higher than that of Yu et al., mainly because we detected many pairs. Related to the reason of our high detection ratio, the substantial difference is that we use all genes transcribed by RNA pol II including widely expressed (housekeeping) genes, while they only use tissue specific genes. Other differences include the evaluation model, source of known interactions (they used DIP [20] and TRANSFAC), and clustering of motifs. For the evaluation model, we evaluated motif pairs with two statistical tests: (1) chi-square test for biases of their relative positions from TSS and (2) K-S test for their distance distribution biases. Yu et al. evaluated motif pairs from two viewpoints: (1) over-representation of motif pairs and (2) distance distribution biases, which we adopted as the second test. Thus, their model does not consider the distance between the motif and TSS. They evaluated a p-value as the product of the two p-values calculated from the above viewpoints. For the thresholds, we used an FDR of 1% for each test, though 1%, 5%, and 10% are conventionally used [21, 22]. Yu et al. set the threshold to be 10^-6.2 after the multiple testing correction. If we calculate the product of the p-values corresponding to the two tests, their maximum value among the detected motif pairs was 5.9 × 10^-5. In addition, if we employ the same threshold of 10^-6.2, 8,379 out of 8,454 pairs were still significant. Thus, our threshold of FDR is stringent enough.

The Class 1 motif pairs are localized near the TSSs, and may represent components of core promoters [1, 24]. These motif pairs mainly consist of LP and SP types of single motifs, which is notable for two reasons. First, in significantly co-occurring motif pairs, NP types of single motifs are overrepresented (Figure 3). Second, in the heat map that shows the significance of correlated co-occurrence (Figure 1B), the overrepresentations of the co-occurrence of the motif pairs cannot be explained as a direct result of the independent occurrence of each motif around the TSS. Thus, it is possible that the TFs that bind to these motif pairs interact with each other.

The possibility that the Class 1 motif pairs are largely related to constitutive or ubiquitous gene expression is indicated by the three following evidences: (1) the EST analysis (Figure 4), (2) the functions of the TFs that binds to the top 20 motifs that are most frequently used in the Class 1 pairs, and (3) a recent work conducted by another group [25]. For example for the functions of the TFs, the transcription factor E2F regulates gene expression in the cell cycle; SP1, MAZ, NFY, KROX, SP3, and CP2 are the TFs associated with ubiquitous transcription; and PAX3 (rank number 6) is a member of the paired box (PAX) family of transcription factors, that has multiple roles during fetal development. As the third evidences, according to a recent study on alternative promoters [25], single promoters that are CpG-rich and do not have alternative promoters are significantly associated with ubiquitous gene expression. It is likely that this category of promoters is mainly composed of our Class 1 motif pairs.

The relative positions of motif pairs are much more flexible in Class 2 than in Class 1. These pairs appear at least 3 kb away from the TSS, yet most previous studies only examined 1 kb of sequence [26]. These Class 2 motif pairs are mainly composed of NP-type motifs, which are likely to cooperate with other motifs to be functional. In S. cerevisiae, Yu et al. [27] found 300 significant TF interactions (as motif pairs), most of which were not constitutively active. This is consistent with our study of human promoters because the motifs related to constitutive gene expression are thought to be linked to Class 1, which is the minority of the detected motif pairs. Moreover, the top 20 most frequent motifs in Class 2 (as shown in Table 1) appear to be involved in various functions. We show some examples of TF pairs, corresponding to detected motif pairs, that are known to function cooperatively. We detected a motif pair consisting of GATA4 and AP1. GATA4 is thought to regulate genes involved in embryogenesis and in myocardial differentiation and function, and GATA4 and AP1 cooperatively regulate transgenic mice overexpressing cardiac calsequestrin [28]. We have another motif pair consisting of GATA4 and TBX5. TBX5 plays a role in heart development and specification of limb identity together with GATA4 [29]. One of our hub motifs MYOD pairs with several motifs known to be involved in muscle and heart gene regulation (MEF2, SRF, SRF, and AP1) [30]. The motif pair NF-kappaB and C/EBP was also detected, and it is known that these two motifs synergistically regulate the mouse serum amyloid A gene expression induced by inflammatory cytokines [31]. Thus, the motif pairs in Class 2 appear to be involved in various specific functions.

We observed that relative positions of motif pairs in Class 2 ranged widely within the 3 kbp regions. It may be due to alternative promoters [13] or the diversity of TSSs [3]. Kimura et al. [13] showed that at least 52% of human genes have alternative promoter sequences separated by >500 bp. Even with a much more strict threshold of 2 kb, 44% (6,485 of 14,628 genes) were regarded to have alternative promoters (see supplementary Figure 1 of [13]). The other possible reason is that these motif pairs can work at any distance as components of enhancers (cis-regulatory modules). It is well accepted that the location of enhancers is rather arbitrary and that both promoters and enhancers share modular architecture [32]. Thus, it seems likely that some promoters contain enhancer-like modules (in addition to the core region) and that these modules contain the Class 2 motif pairs. In any case, the fact that the top 20 motifs in both classes do not overlap (except for a single motif) supports the view that the two classes are distinct in terms of the way they function.

We selected a set of TFs that work as hubs in the TF-TF interaction network in each class. Table 2 summarizes the overrepresented domains in these hubs. It seems that these domains are commonly characterized by an excess of a certain amino acid(s). This tendency is stronger in Class 2-related TFs than Class 1-related TFs. These biased domains, such as the proline-rich domain, may play common important roles in gene regulation. In fact, the SH3 and WW domains recognize proline-rich peptides. The binding of SH3 domains to proline-rich regions causes the formation of a large number of protein complexes [33]. Most of the proteins that interact with SH3 domains contain at least one PxxP motif [34]. It was recently found that 49% of the entire sequence of human TFs contains intrinsically disordered (ID) regions [35] that differ from DNA binding domains. Using techniques such as NMR (Nuclear Magnetic Resonance), it has been revealed that the trans-activation domains are unstructured in unbound TFs and become structured upon binding to their partners [35]. The ID regions tend to consist of alanine, glutamic acid, glycine, lysine, proline, glutamine, arginine, and serine residues [36]. Among these eight amino acids, five are overrepresented in Table 2. Furthermore, Haynes et al. reported that ID regions are common to hub proteins from four eukaryotic interactomes (C. elegans, S. cerevisiae, D. melanogaster, and Homo sapiens) [37], indicating the importance of ID regions in protein interactions. As an unstructured protein domain with repetitive sequence, the C-terminal domain (CTD) of RNA polymerase II is comprised of a variable number of tandem heptapeptide repeats. This feature is important in the efficient capping, splicing, and polyadenylation of mRNA transcripts [38]. These previous studies lead us to a model in which repetitive domains have potential in terms of acquiring multiple functions through protein interaction. Thus, it is likely that the domains commonly observed in hub-TFs (see Table 2) are related to ID regions; this may play an important role in TF-TF interaction in transcription.

With the intention of identifying general features of promoter architecture, the present study examined a general set of genes rather than a set of co-regulated genes, as done in previous studies [9, 10, 26]. Further studies that employ this approach will help to clarify these general features. A possible approach would be the identification of cis-regulatory modules (CRMs) based on the detected motif pairs in this paper, as attempted in [39]. Given that the unit length employed in the present study was 100 bp, a necessary future study would involve a more detailed analysis of the distances between significant motif pairs. We found that most of the motif pair distances are less than 100 bp, with few in the range of 200 bp. This finding is consistent with the observation that most of the cooperative motifs in the TRANSCompel database (a database of composite promoter elements [14], ver. 9.2) are within a distance of 140 bp. Another necessary future task is a more detailed characterization of the remaining significant motif pairs that we tentatively termed Class 3. It seems certain that our visualization method, as well as our statistical tools, will be useful in obtaining further insights into promoter structures.

Our results indicate the occurrence of at least two distinct classes of motif pairs in human promoters, with one class appearing around TSSs and the other found at somewhat arbitrary distances from TSSs. This trend is reminiscent of the well-known difference between promoters and enhancers (CRMs). There has been some discussion that promoters and enhancers are alike in terms of their modular structure as clusters of TFBSs [32]; however, a recent study reported differences in terms of the pattern of histone modifications [40]. Thus, our findings may reflect this kind of difference or merely differences between core promoters and CRM-like modules of promoters; further investigations are necessary to clarify these hypotheses. Our visualization method will also be useful in these types of future studies.