Annotating significant pairs of transcription factor binding sites in regulatory DNA

---

Abstract

In the presented work we search for transcription factor binding sites (BS) by including additional information about typical BS patterns. The new proposed score combines the ordinary profile score based on TRANSFAC-matrices together with a score based on pairs of BS. The latter score positively weights pairs of BS that tend to occur together in many regulatory DNA-sequences, in contrast to a random background model. The empirical BS pair frequencies result from our evaluation of a large dataset of orthologous genes.

Key words: recognition of genes and regulatory elements, transcription factor binding site, score matrix

---

Introduction

Vital processes in cells require a coordinated regulation of gene transcription. Gene transcription is controlled by protein complexes embedding transcription factors which bind to specific sequence patterns in regulatory DNA. These transcription factor binding sites (BS) occur in combinations of core sites and sites that are specific for groups of genes in consideration of order, distance and strand orientation. Therefore investigating the organization of BS in regulatory DNA is essential for understanding transcriptional regulation.

Binding sites are patterns of conserved and non conserved nucleotides, and can be described by TRANSFAC-matrices [1]. Searching BS with TRANSFAC-matrices in regulatory DNA leads to many putative sites covering biological BS and random hits. A current line of research analyzes combinations of BS to improve the discrimination of biological signals from noise. Pairs of BS have been identified in expression data [2, 3, 4] and in tissue data [5]. Combinations of more than two BS have been studied by refining an initial model of a BS-combination [6, 7, 8] and by applying a genetic algorithm [9]. More general approaches applicable to all transcription factors and to whole genomes have been analyzed in [10, 11]. In [10] pairs of BS in a fixed window are scored against a random background model and in [11] combinations of BS are scored according to their distance.

Today many genomes have been completely sequenced which allows to apply the concept of phylogenetic footprinting on orthologous genes of different organisms [12]. This concept identifies DNA with regulatory function as conserved subsequences in noncoding DNA due to the conjecture that for orthologous genes the transcriptional regulation should be conserved. Thus conserved BS in noncoding DNA have a higher significance for biological function than nonconserved BS. This is an advantage in higher organisms where the noncoding DNA can have a length up to 15kB and conserved subsequences form about 10 percent. Combinations of BS in conserved regulatory DNA have been analyzed in [13, 14, 15, 16].

In this work we search for significant pairs of BS in a large sequence set of conserved regulatory DNA of man and mouse. We obtain the set of significant pairs of BS by the following steps: First we label conserved BS in these sequences with TRANSFAC-matrices. Then we search in this set for pairs. We propose a score to quantify which pairs of BS prefer to occur together in comparison with a random background model. As a result we obtain a set of significant pairs of BS which is of general validity for human and mouse genes. Finally we apply our results to annotate significant pairs of BS in regulatory DNA with labelled putative BS. The annotated significant pairs are stronger signals for biological function than the single sites.

Significant pairs of BS in the whole human genome have been analyzed also in [10]. However, our approach differs by the following properties: We search for pairs in conserved DNA, we take into account pairs which are formed by the same factors, we perform the counting of pairs such that pairs are significant for many genes, we use a different definition of the scores for binding sites and pairs of binding sites and we apply our results to annotate significant pairs of BS in regulatory DNA. We describe our approach in detail under "Data and method". Under "Annotation of significant pairs", we apply our approach to muscle-specific genes [13]. At the end, conclusions will be drawn from the results obtained.

---

Data and method

Transcription factor binding sites

The TRANSFAC database [1] contains matrices which describe profiles of known BS of transcription factors. With these matrices one can search for putative BS in regulatory DNA. A quantitative measure for the similarity between a sequence pattern s = (s₁,...,s_n) in the regulatory DNA and a TRANSFAC-profile is the score S_BS defined as

(1)

where p_i(s_i) is the probability for nucleotide s_i at position i based on the TRANSFAC-profile and b(s_i) is the background probability for nucleotide s_i. The b(s_i) have been assumed to be equidistributed so that the score is independent from the total number of sequences and can also be compared with scores in sequences which do not belong to our sequence set ³.

The score S_BS can be transformed to a score

(2)

with zero mean and unit variance for all matrices. For each matrix and have been computed from a random background sequence with an equidistribution for the nucleotides probabilities. This roughly allows to compare scores of different matrices.

Conserved sequences and putative binding sites

Our sequence set is composed of conserved regulatory sequences of man-mouse orthologs from the CORG database [17]. First we removed redundant sequence-pairs, i. e. all pairs where a gene occurred several times. Then we searched with all human and mouse TRANSFAC-matrices (release 6.2) for BS in the resulting 7574 sequence-pairs. We collected candidate sites which score higher than a fixed threshold according to Eq. (2) and which are conserved between man and mouse. The threshold has the value 4. For this threshold we found the majority of a set of biological sites in human muscle genes as shown in Table 1. The number of putative BS increases by more the twofold for a threshold of 3.8.

Pairs of binding sites

Next we search in the collection of candidate sites for pairs of BS. We are interested in significant coocurrences of pairs. We therefore introduce a score as

(3)

with

and

(4)

where #pair_{i, j} is the number of pairs of the BS i and j, m_ij is the pair probability and π_i and π_j are the associated marginals. A positive score S_ij > 0 indicates that the two factors i and j prefer to occur together.

We search the pairs of BS in a window of 10 to 50 base pairs⁴. The minimal distance of 10 bp is necessary to avoid pairs formed by overlapping BS. We have chosen as maximal distance 50 bp due to the fact that biological BS frequently form closely neighbored clusters. We are interested in pairs which regulate many genes. Thus we count the number of genes in which the pair has been found instead of counting the total number of occurrences of a pair. By this way of counting we avoid that a highly scoring pair could result from a high number of pairs in only one gene. Furthermore we take into account that a factor could prefer to occur without a neighbor. These frequencies divided by two are stored in the last column and last row of the matrix m_ij. We distinguish the order of the TFs in a pair, i. e. upstream–TF₁–TF₂–downstream is another pairs as upstream–TF₂–TF₁–downstream. The strand orientation is neglected. The pair frequency m_{i, j} has been annotated only for pairs occurring more often than five times to avoid a high value for S_ij from a low pair frequency and low marginals.

Figure 1 shows a score-matrix which summarizes the scores S_ij for all pairs of transcription factors obtained from our analysis for the following parameters: organism human, conserved BS, threshold for S_BS is equal 4, ordered pairs, window 10 to 50 base pairs. The dark blue background is related to pairs observed less than six times which form 99 percent of all possible pairs. Ninety percent of 580 frequently occurring pairs have a score S_ij > 0. The figure depicts that some transcription factors form pairs with only one or few other factors, like SRF. Other transcription factors form pairs with many different factors, like MAZ, MEF2 and POU3F2. Frequently a pair is formed by twice the same factor. The twin-pairs (SRF,SRF), (NRF1,NRF1) and (CHX10,CHX10) belong to the pairs with the highest scores S_ij > 2. The BS in all pairs are nonoverlapping because their minimal distance is 10 bp. Looking at regulatory sequences annotated with twin-pairs revealed that they partially have three, four or more than ten binding sites of this TF. However, such clusters of BS exists also for a TF where the twin-pair has a low score as in the case of (MYOGENIN,MYOGENIN). The names for the factors we use in our paper are exactly the names of the TRANSFAC-matrices, because in several cases a matrix describes binding site profiles of different factors.

The lines in Figure 2 connect all pairs of transcription factors with S_ij > 1. Ten pairs which are connected by red colored lines are known interacting factors. For S_ij > 0 we found in total 40 pairs of known interacting factors. The majority of factors is connected by black lines showing that our method has identified many new preferred pairs. Several factors form pairs with different partners and in some cases the partners form pairs too. Thus one could also check for preferred triples. Note that in these figures there is still redundancy because some TRANSFAC-matrices have similar profiles which could be removed by clustering of the matrices.

Comparison of pair scores for different parameters

In this section we compare the results for the pair score S_ij for different sets of parameters.

The top-ranking pairs for the parameters: organism human, conserved BS, threshold for S_BS is equal 4, ordered pairs, window 10 to 50 bp which were discussed in the previous subsection are listed in the first column of Table 2.

The top-ranking pairs in the nonconserved noncoding DNA (281 sequence pairs) are different from the top ranking pairs in the conserved DNA as shown in column two and three of Table 2. This is a further evidence that top-ranking pairs in the conserved regulatory DNA have a conserved regulatory function because pairs in nonconserved DNA are less significant for regulatory function. However, top-ranking pairs in nonconserved noncoding DNA could be formed by organism specific factors, like MTF1 for man genes and COUP for mouse genes but TRANSFAC indicates activity of these TFs in both organisms. An exception is the pair (NRF1,NRF1) which belongs to the top-ranking pairs in all columns⁵. Note that two of the nine pairs of (NRF1,NRF1) in the nonconserved regions occur in the same genes of man.

The top-ranking pairs disregarding the BS-order⁶ are shown in the first column of Table 3. One partially recovers the pairs of the first column of table one but also many new pairs. The pair score S_ij is lower than above. Only fourteen percent of 1064 annotated pairs have a positive pair score S_ij. Twin pairs (i,i) are rare in the top-ranking pairs in contrast to above. Neglecting the order of pairs leads to an increase of m_ij for i j and π_i whereas m_ii remains constant. As a consequence S_ii decreases stronger than S_ij.

A TF can have a weaker protein-DNA binding pattern caused by an interaction with other transcription factors. We have sought for such pairs by reducing the BS-threshold to 3.8. Many new top-pairs arise for this lower score-threshold as shown in the second column of Table 3. The best candidate for a weaker DNA binding pattern are the pairs (CREB,ARNT) where the number of pairs increased by a factor 4 from 3 to 12 pairs whereas the marginals increase by a factor 3 and (ARNT,CREB) where the number of pairs increase by a factor 8 from 1 to 8 pairs and the marginals by factors 3 and 4. We have not found any information in pubmed that these two factors can interact.

In addition we have studied whether there are pairs which prefer to occur in a larger distance. Significant pairs in a window 10-200 bp are shown in the third column of Table 3 and contain many new pairs. We have not found a candidate for a preferred larger distance among the top-ten pairs because for all pairs the number of the pairs and the marginals increase with approximately the same rate.

We have compared our sets of significant pairs with the results for significant pairs in the whole human genome [10]. Pairs formed by the same factor which are significant in our results were not taken into account in a definition of a pair in [10]. We have checked whether the significant pairs reported in Table 1 and 2 in [10] occur also in our sets of significant pairs. We recovered only a small number of these pairs in our significant pairs as summarized in Table 4. A reason for these differences could be the different scoring of the BS and the pairs and the different counting of pairs in the two methods.

---

Annotation of significant pairs

First we have annotated pairs of BS with S_ij > 0 in our sequence-set of conserved regulatory DNA. The results are shown on our web-page.

Furthermore our results can be used to annotate pairs of BS in an arbitrary regulatory DNA with labelled putative BS. We show this for muscle specific genes in this section. The threshold for BS has the value 3.8 in the following however the BS-threshold for the significant pairs has the more conservative value 4 according to the results in the section "Pairs of binding sites", see Figure 1 and Table 1. The other parameters are the same as in the section "Pairs of binding sites".

Muscle-specific genes

Wassermann et al. have studied muscle specific conserved regulatory sequences of man and mouse [13]. Their results for BS of muscle specific factors are in good agreement with experimentally confirmed BS and in addition they have predicted new sites.

We have searched for putative BS in these sequences with all TRANSFAC-matrices. Then we have annotated the pairs of BS which are significant in the whole human and mouse conserved DNA shown in Figure 1 in the set of putative BS in the muscle sequences. Our results are summarized in Table 5. The TF AP4 occurs in the original data instead of MYOGENIN in the genes 3, 10, 53 and 60. Frequently the BS of AP4 and MYOGENIN overlap due to similar matrix-profiles. Thus we have noted for these genes the muscle-specific factor in the table. In addition experimentally confirmed BS are marked boldface.

Our method has found known experimental sites and it could predict new sites. We found conserved BS in 25 genes and we could annotate significant pairs in 15 genes. The 22 annotated pairs contain 13 experimentally confirmed BS. Significant pairs of 10 genes are composed of muscle specific factors. The muscle-specific pairs are (SRF,SRF), (MYOGENIN,MYOGENIN), (MEF2,MEF2) and (MEF2,MYOGENIN).

---

Conclusions

In this work we have introduced a refined method to search for BS by calculating a score for pairs of BS. This score positively weights a pair of BS which prefer to occur together in a large sequence set of conserved regulatory DNA of human and mouse in comparison with a random background model. As a result we have obtained a set of significant pairs of BS which contains known interacting pairs, but also many new pairs. It can be used to annotate significant pairs of BS in any arbitrary regulatory DNA with labelled putative BS. This refined searching for BS improves the discrimination between biological sites and random sites.

We have applied our results to annotate significant pairs in muscle-specific genes. The 22 annotated pairs contain 13 experimentally confirmed binding sites. In future work we will improve the scoring of putative BS [18] and will augment the scoring scheme to include information on the position, on strand orientation, and on the length of the intermediate sequence between BS. The challenges are the functional form of the score and the weighting of the individual pieces of information.

Availability: Full list of results at http://www.molgen.mpg.de/~rateit_k/data.html.

---

Acknowledgements

We acknowledge support from DFG in Sonderforschungsbereich 618.