BACKGROUND
gene expression in multicellular eukaryotes varies considerably between tissues and can change dramatically even within the same cell type.
regulation of transcription is one of the fundamental mechanisms for controlling the observed diversity in gene expression  <cit> , and recent studies have underscored the importance of combinatorial regulation by multiple transcription factors  in this regard  <cit> .
progress is also being made towards experimental methods for testing combinatorial regulators on a larger scale in near-physiological conditions  <cit> .
combinatorial regulation can explain, in general, how a relatively small number of tfs can govern gene expression under diverse conditions.
one such example is the regulation of gene expression in immune responses.
pathogen recognition in the vertebrate innate immune system is initially performed by a limited number of pattern-recognition receptors .
the toll-like receptors  are a family of prrs responsible for the recognition of a wide variety of pathogen-associated ligands, such as lipopolysaccharide, viral rna, unmethylated cpg dna and so on.
the recognition of ligands activates signaling pathways leading to the activation of several tfs, such as nf-κb, and irfs.
these tfs are known to induce expression of various genes and evoke pleiotropic immune responses.
although several studies have addressed the importance of combinatorial transcriptional regulation in tlr signaling  <cit> , little is known about which specific combinations of tfs are involved.
various bioinformatics strategies have been used for the prediction of cooperation between tfs.
some studies have used a combination of features, such as co-expression data and protein-protein interactions  <cit> .
other studies utilized a combination of chromatin immunoprecipitation combined with microarray  data and expression data  <cit> .
however, while the integration of heterogeneous experimental data sources is potentially very powerful, in practice such data is too scarce to be of use to a particular tissue of interest, especially for higher eukaryotes like humans and mice.
chip-chip data, in particular, is available for a very limited number of tfs, in a limited number of cell types.
in addition, in the case of de novo predicted regulatory motifs, it might not be known what protein  is binding the motif in question, which restricts the applicability of chip-chip analysis.
for these reasons, a number of studies have focused on identifying combinatorial regulation solely based on predicted transcription factor binding sites .
for example, murakami et al.
used position weight matrices  to predict tfbss on a genomic scale in order to quantify the co-occurrence of regulatory motifs in human promoters  <cit> .
sudarsanam and colleagues used a cumulative hypergeometric distribution to predict regulatory motifs co-occurring on a genome-wide scale in yeast  <cit> .
other studies have described measures for co-occurrence of pairs of motifs as a measure to predict tf synergy  <cit> .
synthetic libraries of promoters have been used to study combinatorial regulation using thermodynamic models  <cit> , and more recently, combinations of oligomers have been used to predict from sequence ep300-bound and crebbp-bound enhancers in three mammalian cell types  <cit> .
a small number of studies have attempted to identify pairs of co-occurring motifs in the promoters of co-expressed genes  <cit> .
however, methods for predicting combinatorial regulation from predicted tfbss are plagued by a number of problems.
these include similarities beween the pwms used to predict tfbss, biases caused by motif over-representation, and difficulty of evaluating the significance of observed co-occurrences using standard statistical tests.
in this study, first we describe a new measure for tfbs pair co-occurrence.
for each pwm pair , we calculate the frequency of motif b in sequences containing one or more a sites, as well as the frequency of motif b in sequences that lack a sites.
we use the ratio of these two frequencies, the frequency ratio , as a measure for co-occurrence.
applying this measure on the tfbss in the genomic set of human and mouse promoters, we observed how co-occurrence tendencies are strikingly different between promoters with high gc content and cpg scores and promoters with low gc content and cpg scores, with the latter having a higher variety in fr values.
we also observed a strong influence of tfbs gc content differences.
based on the above observations, we developed a method for predicting co-regularing pairs of tfs in a set of co-expressed genes.
given the promoter sequences for a set of genes that are co-expressed, we identify motif pairs that co-occur more often than expected.
we use the relative increase in co-occurence in the co-expressed set of genes as an indicator of combinatorial regulation.
our proposed method was designed to overcome the problems associated with previously reported statistics-based measures of co-regulation.
in order to obtain a measure of statistical significance, we compare the observed fr values for pairs of motifs in a set of co-expressed genes with those of sets of genes sampled at random, thus taking into account biases caused by genome-wide co-occurrence tendencies.
we applied our approach to a number of sets of co-expressed mouse genes, and found a number of significantly co-occurring pwms pairs.
importantly, the proposed approach was not biased by tfbs motif over-representation, and could thus detect co-occurrences missed by existing approaches.
for the identified tf pair nf-κb - c/ebpαwe experimentally validated the co-regulation after tlr stimulation in dendritic cells.
since the proposed method does not rely on chip-chip data, it is generally applicable and can complement existing computational methods for discovery of tf co-regulation.
methods
we refer to additional file  <dig> for a workflow of our framework for the detection of co-occurring motifs.
promoter sequences
we used a combination of dbtss data  <cit> , cage data  <cit> , and annotation data from the ucsc genome browser  <cit>  to define transcription start site  positions for both human and mouse genes, as described before  <cit> .
the regions from - <dig> to + <dig> were extracted from the repeat-masked hg <dig> and mm <dig> versions of the human and mouse genome.
for each pair of highly similar sequences  one sequence was removed from our sequence dataset in order to reduce biases caused by duplicated sequences.
position weight matrix dataset
from the transfac  <cit>  and jaspar  <cit>  databases all vertebrate pwms were extracted.
redundancies were removed using tomtom  <cit>  by the following strategy: for each pair of similar pwms  the motif with the lowest information content was removed from our dataset.
pairs were considered in order of increasing tomtom e value.
this resulted in a pwm dataset of  <dig> non-redundant pwms, each representing a group of similar pwms.
for each pwm a score threshold was set in a way that there is about  <dig> hit per  <dig> bps in the mouse promoter sequences.
gc content values of pwms were calculated as the average of the probability of nucleotides c and g over all positions of the pwms.
measure for tfbs co-occurrence: frequency ratio
as a measure of tfbs co-occurrence we introduce the frequency ratio  value.
consider two tfs, tf a and tf b, whose binding preferences are represented by pwm a and pwm b respectively.
given a set of sequences and the predicted sites for both pwms, we calculate the fr, the tendency of sites for tf b to co-occur with those of tf a, as follows.
first, we define seq as the number of sequences containing at least one site for motif a, and n as the number of sites for motif b co-occurring with one or more sites for motif a. from these we calculate frequency, a measure for the number of b sites co-occurring with a sites:
frequencyb|a=nb|a/seqa
likewise, we define frequency as the number of predicted sites for motif b per sequence lacking sites for motif a:
frequencyb|!a=nb|!a/seq!a
where n is the number of b sites in the set of promoters lacking a sites and seq is the number of sequences without a sites.
we calculate the ratio of these two frequency values, fr:
fr=nseqnseq=numberofbsitespersequencehaving ≥1asitenumberofbsitespersequencehavingnoasite
fr is a measure for the tendency of sites of motif b to be present in sequences having at least one a site compared to sequences not having an a site.
in order to limit the bias caused by overlapping sites for pairs of similar motifs, sites for motif b overlapping a sites are not included in n. in the case of homotypic motif pairs , there are obviously no motif b sites in sequences not containing sites for motif a. in this case we define frequency =  <dig>  in this case fr can be interpreted as the average number of a sites in sequences containing at least one a site.
note that the fr measure is not limited to tfbs motifs, but can be used for other sequence motifs and nucleotide oligomers.
micro-array gene expression data
we used micro-array expression data for a large number of human and mouse tissues  <cit> , and for dendritic cells  after stimulation with a number of immune stimuli  <cit>  .
the raw intensity data were processed to calculate robust multi-array average  values.
genes with at least 3-fold differential expression between any pair of samples were picked up.
expression values for each gene were rescaled to mean  <dig> and standard deviation  <dig>  and dimension reduction was performed with principal component analysis.
the gene expression patterns were hierarchically clustered with ward's algorithm and divided into subclusters.
all analyses were performed using r and bioconductor.
definition of cpghigh promoters and cpglow promoters
we based our definition of cpghigh and cpglow promoters on previously reported definitions for cpg islands  <cit> .
for each promoter sequence we defined the ratio of each nucleotide as the number occurrences of that nucleotide divided by the number of non-masked nucleotides in the sequence.
the gc content is then the sum of the c ratio and the g ratio.
the cpg score is the observed ratio of cpg dinucleotides divided by the ratio expected from the ratio of c and g nucleotides.
we defined "cpghigh promoters" as promoters with a gc content ≥ 50% and a cpg score ≥  <dig> .
promoters not meeting these requirements were labeled as "cpglow promoters".
these two conditions gave us  <dig> cpghigh promoters in mouse  and  <dig> cpghigh promoters in human .
detection of over-represented tfbss
as a measure for over-representation of a tfbs motif in a set of sequences we use the over-representation index  as defined by bajic et al.
<cit> .
p-values for ori values were assigned by random sampling sets of sequences of the same size as the set of interest from the genomic set of promoter sequences.
for the estimation of p-value in cpghigh promoters and cpglow promoters the sampling was done from the genomic set of cpghigh promoters and the genomic set of cpglow promoters, respectively.
as p-value threshold for over-representation we used the threshold of  <dig> .
co-occurrence significance in promoters of co-regulated genes
for each co-expressed set of genes containing at least  <dig> genes, we identified over-represented tfbss, as described above.
for each pair of motifs , where a is an over-represented motif and b is any of the  <dig> pwms, we calculated frset in the promoter sequences corresponding to the set of co-expressed genes.
the significance of frset was evaluated using a sampling approach.
in this sampling approach, a set of sequences are sampled from the genomic set of sequences until the same number of b sites and a sites as in the co-expressed cluster is obtained.
from these, a frsampled is calculated, reflecting the genome wide co-occurrence tendencies of the pair of motifs.
this is repeated a large number of times , and the p-value of frset is defined as the ratio of the number of times where frsampled ≥ frset.
generation of artificial and semi-artificial promoter sequences
artificial promoter sequences were constructed by generating sequences of the same length as the actual promoter sequences used in this study, where at each position the nucleotide is decided using a uniform distribution over the alphabet .
semi-artificial sequences were generated by scanning through actual sequences and randomly adding either a g or c to the semi-artificial sequence when a g or c was encountered in the actual sequence; and randomly adding either an a or t when an a or t was encountered.
construction of plasmids, transfection, and luciferase assay
promoter sequences of selected genes were pcr amplified and cloned into pgl3-basic vectors .
sequences from about - <dig> to + <dig> relative to transcription start site were cloned.
kb-tandem reporters were purchased from promega.
complementary dna for tfs was pcr amplified and cloned into pef-bos expression vectors.
the resulting reporter plasmids and tf over-expression plasmids were co-transfected into hek  <dig> cells with prl-tk encoding renilla luciferase  and appropriate signaling molecules with using lipofectamin  <dig> .
at  <dig> hours after transfection, the cells were lysed and subjected to reporter assay according to the manufacture's instruction .
the primers used will be provided upon request.
RESULTS
frequency ratio, a novel measure for co-occurrence of two tfbss: general results and genomic tendencies
as a measure for the co-occurrence of the tfbss for two tfs, tf a and tf b, we propose the frequency ratio, fr .
the fr value is a measure for the tendency of motif b to co-occur with motif a. on a molecular level, it reflects the tendency of tf b to bind the same promoters as tf a, although this does not necessarily imply a direct physical interaction between a and b. cases where fr values are higher than  <dig> reflect a tendency of tf b to bind promoter sequences that also are bound by tf a, while fr values lower than  <dig> reflect a tendency for tf b to bind to promoter sequences not bound by tf a. to avoid biases caused by motif similarities, sites where motifs a and b overlap were discarded before the calculation of the frequency values.
note that fr is not necessarily the same as, or similar to, fr .
using the above definition of fr we calculated the genome-wide fr values for all  <dig>   tfbs motif pairs, in the genome-wide sets of  <dig>  human promoter sequences, and  <dig>  mouse promoter sequences.
a histogram of fr values in the genomic set of mouse promoters is shown in fig.
1a.
although the majority of pwm pairs have fr values close to  <dig> , some pairs have high or low fr values.
similar observations were made for human sequences .
the outliers with large or small fr values indicate the genome-wide tendencies for high or low co-occurrence of sequence motifs, respectively.
these genome-wide tendencies represent reference values to which we will compare the fr values of particular sets of co-expressed genes.
similar sequence motif pairs tend to be co-occurring
next, we analyzed the correlation between fr values and motif-motif similarity.
we used the difference of gc content between pairs of motifs as an indicator of motif similarity.
fig.
1b shows a plot of the difference in gc content between pairs of motifs versus fr for the genomic set of mouse promoter sequences.
this figure clearly shows that motif pairs with a smaller difference in gc content tend to have higher fr values, while motif pairs with different gc content tend to have lower fr values.
a similar tendency was obtained in human promoter sequences .
since we excluded overlapping sites, the tendency of these motifs to co-occur cannot be explained simply by a tendency of sites for similar pwms to overlap with each other.
in semi-artificial promoter sequences, where overall gc content and local gc content fluctuations were identical to those of real promoter sequences, a similar tendency was observed .
on the other hand, in completely artificial sequences with 50% gc content, this tendency was not observed: the vast majority of motif pairs had fr values close to  <dig> .
cpglow promoters have a higher variety of fr values than cpghigh promoters
given the observed influence of gc content on co-occurrence, we decided to separately investigate the tendencies in cpghigh and cpglow promoters.
the genome-wide set of promoter sequences was divided into a set of cpghigh promoters and cpglow promoters , and the fr values in each set were calculated.
figure 1c shows the difference in gc content between pairs of motifs versus their fr values in the genome-wide cpghigh promoter set in mouse.
from this figure we can see that for the cpghigh promoter set, the tendency for motif pairs with a smaller  difference in gc content to have higher  fr values was not observed.
in contrast, for the cpglow promoters  such a tendency was clearly observed.
these trends were also found in human sequences  and semi-artificial promoters sequences .
assuming that the variety in fr values reflects the potential of sequences to encode combinatorial regulation, these results suggest that the regulatory complexity of cpglow promoters is higher than that of cpghigh promoters.
an additional result supporting the notion that cpglow promoters have a higher potential for combinatorial regulation was obtained from the analysis of fr values of  <dig> randomly selected oligomers.
although these oligomers are likely not to be associated with any regulatory motifs, the tendencies of fr values are similar to those we observed for pwm motifs in cpghigh and cpglow sequences .
only few motif pairs have high or low fr values on a genome-wide level
the above observations raise the question to what extent genome-wide fr values are indicative of combinatorial regulation between pairs of tfs on a genome-wide level.
we compared fr values observed in the genome-wide set of promoter sequences with those observed in  <dig> sets of semi-artificial sequences.
to take into account the influence of gc content difference, we divided pairs of pwms into  <dig> bins according to their pairwise gc content difference.
table  <dig> shows all pairs of pwms with high and low tendencies to co-occur on a genome-wide scale .
*  indicates whether the genomic fr is higher or lower than expected.
two important observations can be made.
the first point is that out of the  <dig>  pairs of pwms, only very few have exceptionally high/low fr values .
in mouse promoters  <dig> pairs, and in human promoters  <dig> pairs of pwms had exceptionally high values, while in mouse promoters  <dig> pairs had low values.
in human and mouse cpghigh sequences and mouse cpglow sequences, no pairs with exceptional values were found.
these indicate that the vast majority of fr values in true sequences are within the range of values we can expect to find in semi-artificial sequences lacking any biological meaning.
a second point is that, interestingly, most of the pairs with high values involve the tbp motif , a core promoter motif.
this motif is thought to be typically present in non-cpg island-associated promoters of genes with relatively strictly regulated transcription initiation.
on the other hand, we found the gc-rich and cpg-rich e2f tf motif, which might be indicative of cpg-rich sequences, to have low co-occurrence with a regulatory motif .
in conclusion, the above observations support the hypothesis that the genome-wide variation in fr values, except for those involving a small number of exceptional sequence motifs such as the tata box and gc-box, is mainly a result of sequence variations, and not an indication of genome-wide combinatorial interactions between tfs.
it is important that these genome-wide biases are taken into account by approaches that predict combinatorial regulation in smaller sets of co-regulated genes.
the fr approach allows for detection of co-occurring regulatory motifs in tissue-specific promoter sequences without bias caused by tfbs over-representation
next, we turned our attention to the problem of finding significantly co-occurring motifs in the promoter sequences of sets of co-expressed genes.
we used gene expression data for a large number of mouse tissues  <cit>  to define clusters of co-expressed genes , and applied our approach to each cluster containing at least  <dig> genes .
below, we present some of our findings for a set of  <dig> genes with high expression in mouse liver and kidney.
of these, only  <dig> genes were associated with a cpghigh promoter.
in the following discussion we focused on the cpglow promoters.
the most over-represented motifs in this set of promoters were those for pwms for hnf <dig>  hnf <dig>  and a number of nuclear receptors .
these factors and their importance in liver-specific regulation of transcription have been widely reported  <cit> .
as described in methods, the significance of co-occurrence for each a-b tf pair present in the set of sequences was estimated by sampling sequences from the genomic set of sequences until the same number of a and b sites were obtained, and the p-value of frset was computed.
in the kidney/liver-specific set we found  <dig> co-occurrences with p-values <  <dig> .
in one case, the pwm for foxp <dig> co-occurred significantly with hnf <dig> sites.
foxp <dig> plays a role in the development of various organs, including liver.
in this set of promoters,  <dig> out of  <dig> non-overlapping foxp <dig> sites co-occurred with hnf <dig> sites, yielding a fr value of  <dig>  .
given the genomic fr value , we would expect only about  <dig>  foxp <dig> sites to co-occur with hnf <dig> sites on average .
moreover, the distance between foxp <dig> sites and hnf <dig> sites was biased towards proximal positioning: of the  <dig> foxp1-hnf <dig> site pairs,  <dig> were separated by less than  <dig> bps .
in addition, visual inspection of the site pairs revealed a preference of the foxp <dig> site to be upstream of the hnf <dig> site.
to illustrate the difference between our approach and approaches based on statistical tests, we calculated co-occurrence p-values using the method of yu et al.
<cit> , and using the method of sudarsanam et al.
<cit> .
the approach by yu et al.
evaluates co-occurrences using two p-values, one for co-occurrences, pocc, and one for the bias in distances between pairs of sites, pd.
here we focused on pocc, the probability of observing an equal or greater number of co-occurrences, calculated based on the number of sequences in the co-regulated set versus the size of the genome-wide set, the number of co-occurrences between two motifs in the genome-wide set, and the number of co-occurrences in the co-expressed set.
the approach by sudarsanam et al.
uses a cumulative hypergeometric model to evaluate the significance of the observed number of co-occurrences for a motif pair, by comparing it to the distribution of expected co-occurrences given the number of occurrences of the individual motifs.
we applied our fr approach, the pocc approach, and the sudarsanam approach on all sets of co-expressed genes, and compared the results in terms of the over-representation of co-occurring motifs.
fig.
<dig> shows that the distribution of ori p-values for all  <dig> pwms co-occurring significantly with an over-represented motif is similar to that of all pwms, confirming that the fr approach is not biased by motif over-representation.
indeed, the majority of predicted co-occurring motifs are not over-represented.
in contrast, the distribution of ori p-values of predicted co-occurring motifs in the top  <dig> pairs as predicted by pocc, showed a strong bias towards lower ori p-values, indicating that this method is strongly biased by motif over-representation.
the fact that with increasing motif over-representation the expected number of co-occurrences modeled by the hypergeometric distribution also increases, makes the approach described by sudarsanam et al.
<cit>  relatively robust against the bias caused by motif over-representation, but less so than the fr measure.
however, this method does not use a reference set of sequences during the evaluation of significance, making it the most easily affected of these three approaches by pwm-to-pwm similarities .
a relatively high number of co-occurring pairs predicted by the approach by sudarsanam et al., have similar gc content levels, and pairs of motifs with large differences in gc content are relatively rarely predicted to be co-occurring .
as an illustration, for the set of promoters of liver- and kidney-specific genes in mouse, the top co-occurrences in terms of pocc were strongly dominated by pwm pairs containing hnf <dig> and hnf <dig>  which were both strongly over-represented in this cluster.
in the top  <dig> motif pairs,  <dig> involved hnf <dig>  which was found to have significant pocc values with most other over-represented motifs, such as those for hnf <dig> and ikaros.
the pair hnf <dig> - hnf <dig> had the lowest pocc value .
however frset was only moderately higher than frgenomic .
indeed, only  <dig> out of  <dig>  hnf <dig> sites co-occurred with hnf <dig> sites, which were present in  <dig> out of  <dig>  sequences in this cluster.
even though both motifs were over-represented in this cluster, they did not have a strong tendency to be present in the same sequences.
the measure described by sudarsanam et al.
too, predicted a number of significant co-occurrences involving hnf <dig>  strikingly, the top  <dig> motifs predicted to co-occur with hnf <dig> motifs have similarly low gc content value as the hnf <dig> motif , while this is not the case for the  <dig> motifs with most significant fr values with regard to hnf <dig> .
collectively, these results indicate that pocc is more related to co-over-representation than to co-occurrence, and that the measure proposed by sudarsanam et al.
is relatively sensitive to gc content similarities.
our proposed approach, on the other hand, can find non-over-represented motifs, which are likely to be missed by traditional approaches, and is less influenced by pwm-to-pwm gc content similarity.
identification of tfs co-occurring with nf-κb or irf in gene sets having specific expression patterns on tlr stimulation
several tfs such as nf-κb, irf, and ap- <dig> are known to play roles in gene expression evoked by tlr signaling.
however, little is known about involvement of other tfs and how these tfs orchestrate variety of gene expression patterns.
to elucidate the diverse gene expression patterns generated by combinations of various tfs, we next applied our method for the analysis on tlr signaling.
we used a large scale microarray dataset on gene expression in bone marrow-derived dcs after stimulation with various tlr stimuli  <cit> .
the gene expression values were calculated and hierarchically clustered into  <dig> clusters .
frequency ratio analysis revealed that a wide variety of tfs had significant co-occurrence with each other .
firstly, we looked into tfs co-occurring with the binding motifs of nf-κb or irf, since ori analysis showed that these two motifs are significantly over-represented in the clusters, suggesting their key roles in tlr signaling, consistent with former studies .
several tfs were found to have significant co-occurrence with nf-κb or irf .
nf-κb and irf motifs were found to co-occur significantly with themselves  in clusters  <dig>   <dig>  and  <dig>  clusters  <dig> and  <dig> contained genes whose expression levels peaked at about  <dig>  to  <dig> hours after stimulation, suggesting multiple nf-κb sites in promoters of those genes enable rapid induction .
genes in cluster  <dig> were specific for poly  stimulation, and their expression peaked at around  <dig> hours after stimulation .
in addition to various irf family tfs, the pwm m <dig> represents stat <dig>  implying that these genes may be induced secondarily by type i interferon.
other well-known motifs found were c-fos , presumably representing the ap- <dig> motif.
these collectively suggest that our method could identify known tf pairs involved in tlr signaling.
we found that in cluster  <dig>  ccaat enhancer binding protein alpha  had significant co-occurrence with nf-κb .
c/ebp family transcription factors are reported to be involved in tlr signaling-induced gene expression such as cytokine gene expression  <cit> .
importantly, while pocc is relatively low for the pair nf-κb - c/ebpα, up to  <dig> out of  <dig> pwms have a pocc <  <dig>  for co-occurrence with the over-represented nf-κb motif in cluster  <dig> , making pocc not useful for predicting co-occurrences in this case.
a gene in the cluster, nos <dig>  has c/ebp motifs and nf-κb motifs in its promoter , thus the promoter activity is expected to be controlled by c/ebp and/or nf-κb.
to check this, we cloned the nos <dig> promoter, and its activity upon over-expression of tfs or signaling molecules myd <dig> simultaneously with c/ebpα was checked by luciferase assay.
the activity of the nos <dig> promoter was up-regulated only by over-expression of c/ebpα, rela , or myd <dig>  , indicating these tfs positively regulate the nos <dig> promoter.
moreover, when c/ebpα was over-expressed simultaneously with rela or myd <dig>  luciferase activity increased compared to that on rela or myd <dig> over-expression alone.
this result indicated that c/ebpα controls the expression of the nos <dig> gene, and also suggested that it controls the expression of genes other than nos <dig> in clusters  <dig>
synergistic activation of tlr-regulated promoters by nf-κb and c/ebpα
the above results prompted us to check if a broader array of promoters is regulated by c/ebpα.
we tested whether c/ebpα controls nf-κb-regulated promoters or not.
since four pwms in our pwm set represent c/ebp tfs , we picked up promoters having predicted tfbss for nf-κb and one of the four motifs .
we also added as a positive control arg <dig>  which has a predicted nf-κb site and has been reported to be a target of c/ebpβ  <cit> .
five promoters out of  <dig> tested were up-regulated synergistically by rela and/or myd <dig> and c/ebpα , whereas activation of one promoter was suppressed by c/ebpα over-expression .
in contrast, the tandem-kb luciferase reporter was suppressed by c/ebpα over-expression , indicating a specific activation of promoters by c/ebpα.
these results suggested that c/ebpα synergistically and specifically up-regulates the activity of a set of promoters regulated by nf-κb.
taken together, our method successfully identified a pair of transcription factors involved in the immune response.
CONCLUSIONS
in this study, we introduced a new measure for regulatory motif co-occurrence, and investigated genome-wide co-occurrence tendencies between pairs of regulatory motifs.
our initial results show that some motif pairs have a strong tendency to co-occur, while other pairs have a strong tendency to avoid co-occurrence.
however, further investigation showed that these tendencies reflect gc content fluctuations in promoter sequences, rather than a genome-wide level of combinatorial regulation: semi-artificial sequences in which gc content fluctuations were identical to real sequences showed very similar trends.
on the other hand, this trend was lost in completely artificial sequences.
also, we found that only few pairs of regulatory motifs had exceptionally high or low fr values in the genomic set of sequences as compared to semi-artificial sequences.
pairs that did show exceptionally high fr values often involved the tata-box motif, which might reflect the tendency of strictly regulated non-cpg island-associated promoters to contain a tata-box more frequently than other promoter sequences.
in addition, we showed that the tendencies were completely different between cpghigh promoters and cpglow promoters.
our measure for regulatory motif co-occurrences showed a relatively limited range in cpghigh promoters compared to cpglow promoters.
similar observations were made for a set of randomly selected oligonucleotide motifs.
these observations reflect a fundamental difference between these two types of promoters.
cpg islands have been reported to be associated with ubiquitously expressed genes and housekeeping genes, while genes not associated with cpg islands tend to be tissue-specific or condition-specific genes  <cit> .
on the promoter sequence level too, there are considerable differences: while the promoters of cpg island-associated genes tend to lack typical core promoter elements and tissue-specific tfbss, the promoters of genes not associated with cpg islands tend to contain tata boxes or other core promoter elements and tfbss allowing their precise regulation of expression  <cit> .
the apparent lack of high and low fr values in the genome-wide set of cpghigh promoters might reflect a relatively low need for complex combinatorial regulation, compared to cpglow promoters.
on the other hand, for tissue- or condition-specific genes combinatorial regulation might be necessary to ensure spatio-temporal specificity, reflected in the larger range of fr values observed in cpglow promoters.
keeping the above observations in mind, in the proposed method, for the set of promoter sequences of interest, the significance of co-occurring pairs was estimated using a random sampling procedure.
this approach thus takes into account the genomic tendency of motif pairs of similar structure to appear in the same promoter sequence.
furthermore, we considered cpghigh promoters and cpglow promoters as separate cases.
a recent study on tfbs analysis has led to a similar recommendation  <cit> .
moreover, we excluded overlapping pairs of sites to avoid bias caused by similarity between motifs.
applying our method to a large number of tissue-specific sets of mouse promoters, we could predict a large number of pairs of significantly co-occurring tfbs pairs.
one example is the pair hnf <dig> - foxp <dig>  for which we found binding sites to be significantly co-occurring in the promoters of genes with specific expression in liver and kidney.
moreover, the tfbss of this pair of tfs showed a tendency to be located proximally to each other, with the foxp <dig> tfbss located upstream of the hnf <dig> tfbss.
importantly, our approach demonstrated improved robustness against biases caused by strongly over-represented motifs in comparison to a previously reported statistics-based method.
indeed, the majority of significant interactions we found involved motifs that were not over-represented .
this was also the case for c/ebpα binding sites in cluster  <dig> from the dc expression data.
such motifs would thus not be detected by standard over-representation analysis.
for one of the significantly co-occurring tf pairs involved in tlr signaling we could verify the predicted combinatorial regulation.
we found that c/ebpα co-regulates a set of promoters with nf-κb.
co-regulation by nf-κb and c/ebp has previously been reported.
lcn <dig> and arg <dig> were reported as targets of c/ebpβ  <cit> , confirming that the frequency ratio could predict biologically meaningful tf pairs.
moreover, c/ebpα itself has recently been found to control a number of rela-dependent inflammatory promoters, and nf-κb activation synergistically with pu.
<dig>  <cit> , further supporting our findings on the significance of the nf-κb-c/ebp pair.
it would be interesting to check the involvement of other pairs identified as in table  <dig> in tlr-induced gene expression patterns.
there is some room for possible improvements of our approach.
first of all, epigenetic factors that might be responsible for tissue- or condition-specific expression should be taken into account.
at present, chromatin remodeling data is still limited to a small number of cell types, which makes it difficult to incorporate in our approach.
undoubtedly, as the amount of available data increases, there will be a need to incorporate it, resulting in approaches combining both tfbs information and epigenetic information.
importantly, while we have limited our analysis to promoter sequences here, the role of distal enhancers in the regulation of transcription is generally accepted.
as epigenetic data for various cell types increases, we will become able to apply our method not only to promoter regions but also to enhancers, and investigate potential differences in combinatorial regulation occurring in promoter and enhancers.
secondly, since our approach relies on tfbs prediction, which is still known to have a low specificity, further developments in the prediction of tfbss, and additional genome-wide binding data is likely to improve our method as well as other tfbs prediction-based methods.
in conclusion, our fr approach circumvents biases which former methodology suffers from, and we could identify some meaningful co-occurring tfbs pairs, one of which was experimentally supported.
we believe this approach can help us detect combinatorial interactions between tfs in the regulation of transcription, and we also believe that this sets a basis for future developments in computational identification of combinatorial gene regulation.
an online application of our method, which we call regulatory motif combination detector , is available at our website  <cit> .
list of abbreviations used
c/ebp: ccaat enhancer binding protein; chip: chromatin immunoprecipitation; dc: dendritic cell; fr: frequency ratio; ori: over-representation index; ppr: pattern-recognition receptor; pwm: position weight matrix; rma: robust multi-array average; tf: transcription factor; tfbs: transcription factor binding site; tlr: toll-like receptor; tss: transcription start site.
competing interests
the authors declare that they have no competing interests.
authors' contributions
av and yk conceived of the study.
av performed the bioinformatics analysis, and prepared the manuscript.
yk carried out experimental validations and prepared the manuscript.
dms and sa supervised the project and helped with discussion and interpretation of results and with drafting the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
figure s <dig> -  workflow of our framework for the detection of co-occurring motifs.
the analysis of genome-wide tendencies starts with a set of tfbss, predicted in promoter sequences and a set of pwms.
for each pair of motifs, fr values are calculated, and used for further analysis of genome-wide tendencies.
the analysis of co-occurrences in sets of co-regulated genes similarly starts with the prediction of tfbss.
using these, significantly over-represented tfbss are detected, and for each motif the tendency to co-occur with each of the over-represented motifs is analysed.
the significance of the co-occurrences is evaluated using a random sampling approach, sampling sequences from the genomic set of promoters.
click here for file
additional file 2
supporting text -  on the asymmetry of the frequency ration measure.
click here for file
additional file 3
figure s <dig> -  genome-wide tendencies of frequency ratios in human promoter sequences.
histogram of fr values for all pwm pairs in the genomic set of human promoter sequences.
plots of gc content differences as measure of pwm-to-pwm dissimilarity  versus fr values , for all promoters , cpghigh promoters , and cpglow promoters .
click here for file
additional file 4
figure s <dig> -  tendencies of frequency ratio in semi-artificial and completely artificial sequences.
plot of gc content differences as measure of pwm-to-pwm dissimilarity  versus fr values  in semi-artificial sequences , and completely artificial sequences , semi-artificial cpghigh sequences , and semi-artificial cpglow sequences .
click here for file
additional file 5
figure s <dig> -  genome-wide tendencies of frequency ratios for  <dig> randomly selected 7-mers in human and mouse promoter sequences.
plots of gc content differences  versus fr values  are shown for all human promoters , all mouse promoters , human cpghigh promoters , mouse cpghigh promoters , human cpglow promoters , and mouse cpglow promoters .
click here for file
additional file 6
figure s <dig> -  heatmap representation of the average expression values for each of the  <dig> clusters obtained from the gnf geneatlas mouse data.
click here for file
additional file 7
table s <dig> -  summary of main tissues for the  <dig> clusters obtained from the gnf geneatlas data.
click here for file
additional file 8
table s <dig> -  summary of over-represented pwm motifs in tissue-specific sets of mouse promoters
click here for file
additional file 9
figure s <dig> -  histogram of the pwm-to-pwm gc content differences of co-occurring motifs predicted by three approaches.
co-occurrences predicted by the fr measure are least affected by pwm-to-pwm gc content differences.
the distribution of gc content differences of predicted co-occurring pairs of pwms is shown 1) for the  <dig> pwms found to be significantly co-occurring with an over-represented motif according to fr values , 2) for the pwms found to be co-occurring with an over-represented motif according to pocc , and 3) for the pwms found to be co-occurring with an over-represented motif according to the approach of sudarsanam et al.
.
for the latter two approaches the  <dig> pairs with the most significant co-occurrence were used.
click here for file
additional file 10
figure s <dig> -  heatmap representation of clusters of tlr-stimulated dc gene expression data referred to in the main text.
click here for file
additional file 11
table s <dig> -  summary for the co-occurrences in tissue-specific sets of mouse promoters .
click here for file
acknowledgements
the authors would like to thank all members of the laboratory of systems immunology and the laboratory of host defense for helpful discussions and advice; karlou mar s. amada for help constructing the online tool; dr k. nakai for helpful suggestions; e. kurumatani, a. yoshimura, and y. kimura for technical assistance; m. ogawa, e. kamada and m. kageyama for secretarial assistance.
computational time was partly provided by the super computer system at the human genome center, institute of medical sciences, the university of tokyo.
funding: this research is supported by the japan society for the promotion of science  through the "funding program for world-leading innovative r&d on science and technology ," initiated by the council for science and technology policy , and by a kakenhi grant-in-aid for scientific research  from the japan society for the promotion of science.
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2012: eleventh international conference on bioinformatics : computational biology.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s <dig>

 gene expression in multicellular eukaryotes varies considerably between tissues and can change dramatically even within the same cell type. regulation of transcription is one of the fundamental mechanisms for controlling the observed diversity in gene expression, and recent studies have underscored the importance of combinatorial regulation by multiple transcription factors in this regard. progress is also being made towards experimental methods for testing combinatorial regulators on a larger scale in near-physiological conditions. combinatorial regulation can explain, in general, how a relatively small number of tfs can govern gene expression under diverse conditions. one such example is the regulation of gene expression in immune responses. pathogen recognition in the vertebrate innate immune system is initially performed by a limited number of pattern-recognition receptors. the toll-like receptors are a family of prrs responsible for the recognition of a wide variety of pathogen-associated ligands, such as lipopolysaccharide, viral rna, unmethylated cpg dna and so on. the recognition of ligands activates signaling pathways leading to the activation of several tfs, such as nf-κb, and irfs. these tfs are known to induce expression of various genes and evoke pleiotropic immune responses. although several studies have addressed the importance of combinatorial transcriptional regulation in tlr signaling, little is known about which specific combinations of tfs are involved. various bioinformatics strategies have been used for the prediction of cooperation between tfs. some studies have used a combination of features, such as co-expression data and protein-protein interactions. other studies utilized a combination of chromatin immunoprecipitation combined with microarray data and expression data. however, while the integration of heterogeneous experimental data sources is potentially very powerful, in practice such data is too scarce to be of use to a particular tissue of interest, especially for higher eukaryotes like humans and mice. chip-chip data, in particular, is available for a very limited number of tfs, in a limited number of cell types. in addition, in the case of de novo predicted regulatory motifs, it might not be known what protein is binding the motif in question, which restricts the applicability of chip-chip analysis. for these reasons, a number of studies have focused on identifying combinatorial regulation solely based on predicted transcription factor binding sites. for example, murakami et al. used position weight matrices to predict tfbss on a genomic scale in order to quantify the co-occurrence of regulatory motifs in human promoters. sudarsanam and colleagues used a cumulative hypergeometric distribution to predict regulatory motifs co-occurring on a genome-wide scale in yeast. other studies have described measures for co-occurrence of pairs of motifs as a measure to predict tf synergy. synthetic libraries of promoters have been used to study combinatorial regulation using thermodynamic models, and more recently, combinations of oligomers have been used to predict from sequence ep300-bound and crebbp-bound enhancers in three mammalian cell types. a small number of studies have attempted to identify pairs of co-occurring motifs in the promoters of co-expressed genes. however, methods for predicting combinatorial regulation from predicted tfbss are plagued by a number of problems. these include similarities beween the pwms used to predict tfbss, biases caused by motif over-representation, and difficulty of evaluating the significance of observed co-occurrences using standard statistical tests. in this study, first we describe a new measure for tfbs pair co-occurrence. for each pwm pair, we calculate the frequency of motif b in sequences containing one or more a sites, as well as the frequency of motif b in sequences that lack a sites. we use the ratio of these two frequencies, the frequency ratio, as a measure for co-occurrence. applying this measure on the tfbss in the genomic set of human and mouse promoters, we observed how co-occurrence tendencies are strikingly different between promoters with high gc content and cpg scores and promoters with low gc content and cpg scores, with the latter having a higher variety in fr values. we also observed a strong influence of tfbs gc content differences. based on the above observations, we developed a method for predicting co-regularing pairs of tfs in a set of co-expressed genes. given the promoter sequences for a set of genes that are co-expressed, we identify motif pairs that co-occur more often than expected. we use the relative increase in co-occurence in the co-expressed set of genes as an indicator of combinatorial regulation. our proposed method was designed to overcome the problems associated with previously reported statistics-based measures of co-regulation. in order to obtain a measure of statistical significance, we compare the observed fr values for pairs of motifs in a set of co-expressed genes with those of sets of genes sampled at random, thus taking into account biases caused by genome-wide co-occurrence tendencies. we applied our approach to a number of sets of co-expressed mouse genes, and found a number of significantly co-occurring pwms pairs. importantly, the proposed approach was not biased by tfbs motif over-representation, and could thus detect co-occurrences missed by existing approaches. for the identified tf pair nf-κb - c/ebpαwe experimentally validated the co-regulation after tlr stimulation in dendritic cells. since the proposed method does not rely on chip-chip data, it is generally applicable and can complement existing computational methods for discovery of tf co-regulation. methods we refer to additional file for a workflow of our framework for the detection of co-occurring motifs. promoter sequences we used a combination of dbtss data, cage data, and annotation data from the ucsc genome browser to define transcription start site positions for both human and mouse genes, as described before. the regions from - to + were extracted from the repeat-masked hg and mm versions of the human and mouse genome. for each pair of highly similar sequences one sequence was removed from our sequence dataset in order to reduce biases caused by duplicated sequences. position weight matrix dataset from the transfac and jaspar databases all vertebrate pwms were extracted. redundancies were removed using tomtom by the following strategy: for each pair of similar pwms the motif with the lowest information content was removed from our dataset. pairs were considered in order of increasing tomtom e value. this resulted in a pwm dataset of non-redundant pwms, each representing a group of similar pwms. for each pwm a score threshold was set in a way that there is about hit per bps in the mouse promoter sequences. gc content values of pwms were calculated as the average of the probability of nucleotides c and g over all positions of the pwms. measure for tfbs co-occurrence: frequency ratio as a measure of tfbs co-occurrence we introduce the frequency ratio value. consider two tfs, tf a and tf b, whose binding preferences are represented by pwm a and pwm b respectively. given a set of sequences and the predicted sites for both pwms, we calculate the fr, the tendency of sites for tf b to co-occur with those of tf a, as follows. first, we define seq as the number of sequences containing at least one site for motif a, and n as the number of sites for motif b co-occurring with one or more sites for motif a. from these we calculate frequency, a measure for the number of b sites co-occurring with a sites: frequencyb|a=nb|a/seqa likewise, we define frequency as the number of predicted sites for motif b per sequence lacking sites for motif a: frequencyb|!a=nb|!a/seq!a where n is the number of b sites in the set of promoters lacking a sites and seq is the number of sequences without a sites. we calculate the ratio of these two frequency values, fr: fr=nseqnseq=numberofbsitespersequencehaving ≥1asitenumberofbsitespersequencehavingnoasite fr is a measure for the tendency of sites of motif b to be present in sequences having at least one a site compared to sequences not having an a site. in order to limit the bias caused by overlapping sites for pairs of similar motifs, sites for motif b overlapping a sites are not included in n. in the case of homotypic motif pairs, there are obviously no motif b sites in sequences not containing sites for motif a. in this case we define frequency = in this case fr can be interpreted as the average number of a sites in sequences containing at least one a site. note that the fr measure is not limited to tfbs motifs, but can be used for other sequence motifs and nucleotide oligomers. micro-array gene expression data we used micro-array expression data for a large number of human and mouse tissues, and for dendritic cells after stimulation with a number of immune stimuli. the raw intensity data were processed to calculate robust multi-array average values. genes with at least 3-fold differential expression between any pair of samples were picked up. expression values for each gene were rescaled to mean and standard deviation and dimension reduction was performed with principal component analysis. the gene expression patterns were hierarchically clustered with ward's algorithm and divided into subclusters. all analyses were performed using r and bioconductor. definition of cpghigh promoters and cpglow promoters we based our definition of cpghigh and cpglow promoters on previously reported definitions for cpg islands. for each promoter sequence we defined the ratio of each nucleotide as the number occurrences of that nucleotide divided by the number of non-masked nucleotides in the sequence. the gc content is then the sum of the c ratio and the g ratio. the cpg score is the observed ratio of cpg dinucleotides divided by the ratio expected from the ratio of c and g nucleotides. we defined "cpghigh promoters" as promoters with a gc content ≥ 50% and a cpg score ≥. promoters not meeting these requirements were labeled as "cpglow promoters". these two conditions gave us cpghigh promoters in mouse and cpghigh promoters in human. detection of over-represented tfbss as a measure for over-representation of a tfbs motif in a set of sequences we use the over-representation index as defined by bajic et al.. p-values for ori values were assigned by random sampling sets of sequences of the same size as the set of interest from the genomic set of promoter sequences. for the estimation of p-value in cpghigh promoters and cpglow promoters the sampling was done from the genomic set of cpghigh promoters and the genomic set of cpglow promoters, respectively. as p-value threshold for over-representation we used the threshold of. co-occurrence significance in promoters of co-regulated genes for each co-expressed set of genes containing at least genes, we identified over-represented tfbss, as described above. for each pair of motifs, where a is an over-represented motif and b is any of the pwms, we calculated frset in the promoter sequences corresponding to the set of co-expressed genes. the significance of frset was evaluated using a sampling approach. in this sampling approach, a set of sequences are sampled from the genomic set of sequences until the same number of b sites and a sites as in the co-expressed cluster is obtained. from these, a frsampled is calculated, reflecting the genome wide co-occurrence tendencies of the pair of motifs. this is repeated a large number of times, and the p-value of frset is defined as the ratio of the number of times where frsampled ≥ frset. generation of artificial and semi-artificial promoter sequences artificial promoter sequences were constructed by generating sequences of the same length as the actual promoter sequences used in this study, where at each position the nucleotide is decided using a uniform distribution over the alphabet. semi-artificial sequences were generated by scanning through actual sequences and randomly adding either a g or c to the semi-artificial sequence when a g or c was encountered in the actual sequence; and randomly adding either an a or t when an a or t was encountered. construction of plasmids, transfection, and luciferase assay promoter sequences of selected genes were pcr amplified and cloned into pgl3-basic vectors. sequences from about - to + relative to transcription start site were cloned. kb-tandem reporters were purchased from promega. complementary dna for tfs was pcr amplified and cloned into pef-bos expression vectors. the resulting reporter plasmids and tf over-expression plasmids were co-transfected into hek cells with prl-tk encoding renilla luciferase and appropriate signaling molecules with using lipofectamin. at hours after transfection, the cells were lysed and subjected to reporter assay according to the manufacture's instruction. the primers used will be provided upon request.  frequency ratio, a novel measure for co-occurrence of two tfbss: general results and genomic tendencies as a measure for the co-occurrence of the tfbss for two tfs, tf a and tf b, we propose the frequency ratio, fr. the fr value is a measure for the tendency of motif b to co-occur with motif a. on a molecular level, it reflects the tendency of tf b to bind the same promoters as tf a, although this does not necessarily imply a direct physical interaction between a and b. cases where fr values are higher than reflect a tendency of tf b to bind promoter sequences that also are bound by tf a, while fr values lower than reflect a tendency for tf b to bind to promoter sequences not bound by tf a. to avoid biases caused by motif similarities, sites where motifs a and b overlap were discarded before the calculation of the frequency values. note that fr is not necessarily the same as, or similar to, fr. using the above definition of fr we calculated the genome-wide fr values for all tfbs motif pairs, in the genome-wide sets of human promoter sequences, and mouse promoter sequences. a histogram of fr values in the genomic set of mouse promoters is shown in fig. 1a. although the majority of pwm pairs have fr values close to, some pairs have high or low fr values. similar observations were made for human sequences. the outliers with large or small fr values indicate the genome-wide tendencies for high or low co-occurrence of sequence motifs, respectively. these genome-wide tendencies represent reference values to which we will compare the fr values of particular sets of co-expressed genes. similar sequence motif pairs tend to be co-occurring next, we analyzed the correlation between fr values and motif-motif similarity. we used the difference of gc content between pairs of motifs as an indicator of motif similarity. fig. 1b shows a plot of the difference in gc content between pairs of motifs versus fr for the genomic set of mouse promoter sequences. this figure clearly shows that motif pairs with a smaller difference in gc content tend to have higher fr values, while motif pairs with different gc content tend to have lower fr values. a similar tendency was obtained in human promoter sequences. since we excluded overlapping sites, the tendency of these motifs to co-occur cannot be explained simply by a tendency of sites for similar pwms to overlap with each other. in semi-artificial promoter sequences, where overall gc content and local gc content fluctuations were identical to those of real promoter sequences, a similar tendency was observed. on the other hand, in completely artificial sequences with 50% gc content, this tendency was not observed: the vast majority of motif pairs had fr values close to. cpglow promoters have a higher variety of fr values than cpghigh promoters given the observed influence of gc content on co-occurrence, we decided to separately investigate the tendencies in cpghigh and cpglow promoters. the genome-wide set of promoter sequences was divided into a set of cpghigh promoters and cpglow promoters, and the fr values in each set were calculated. figure 1c shows the difference in gc content between pairs of motifs versus their fr values in the genome-wide cpghigh promoter set in mouse. from this figure we can see that for the cpghigh promoter set, the tendency for motif pairs with a smaller difference in gc content to have higher fr values was not observed. in contrast, for the cpglow promoters such a tendency was clearly observed. these trends were also found in human sequences and semi-artificial promoters sequences. assuming that the variety in fr values reflects the potential of sequences to encode combinatorial regulation, these results suggest that the regulatory complexity of cpglow promoters is higher than that of cpghigh promoters. an additional result supporting the notion that cpglow promoters have a higher potential for combinatorial regulation was obtained from the analysis of fr values of randomly selected oligomers. although these oligomers are likely not to be associated with any regulatory motifs, the tendencies of fr values are similar to those we observed for pwm motifs in cpghigh and cpglow sequences. only few motif pairs have high or low fr values on a genome-wide level the above observations raise the question to what extent genome-wide fr values are indicative of combinatorial regulation between pairs of tfs on a genome-wide level. we compared fr values observed in the genome-wide set of promoter sequences with those observed in sets of semi-artificial sequences. to take into account the influence of gc content difference, we divided pairs of pwms into bins according to their pairwise gc content difference. table shows all pairs of pwms with high and low tendencies to co-occur on a genome-wide scale. * indicates whether the genomic fr is higher or lower than expected. two important observations can be made. the first point is that out of the pairs of pwms, only very few have exceptionally high/low fr values. in mouse promoters pairs, and in human promoters pairs of pwms had exceptionally high values, while in mouse promoters pairs had low values. in human and mouse cpghigh sequences and mouse cpglow sequences, no pairs with exceptional values were found. these indicate that the vast majority of fr values in true sequences are within the range of values we can expect to find in semi-artificial sequences lacking any biological meaning. a second point is that, interestingly, most of the pairs with high values involve the tbp motif, a core promoter motif. this motif is thought to be typically present in non-cpg island-associated promoters of genes with relatively strictly regulated transcription initiation. on the other hand, we found the gc-rich and cpg-rich e2f tf motif, which might be indicative of cpg-rich sequences, to have low co-occurrence with a regulatory motif. in conclusion, the above observations support the hypothesis that the genome-wide variation in fr values, except for those involving a small number of exceptional sequence motifs such as the tata box and gc-box, is mainly a result of sequence variations, and not an indication of genome-wide combinatorial interactions between tfs. it is important that these genome-wide biases are taken into account by approaches that predict combinatorial regulation in smaller sets of co-regulated genes. the fr approach allows for detection of co-occurring regulatory motifs in tissue-specific promoter sequences without bias caused by tfbs over-representation next, we turned our attention to the problem of finding significantly co-occurring motifs in the promoter sequences of sets of co-expressed genes. we used gene expression data for a large number of mouse tissues to define clusters of co-expressed genes, and applied our approach to each cluster containing at least genes. below, we present some of our findings for a set of genes with high expression in mouse liver and kidney. of these, only genes were associated with a cpghigh promoter. in the following discussion we focused on the cpglow promoters. the most over-represented motifs in this set of promoters were those for pwms for hnf hnf and a number of nuclear receptors. these factors and their importance in liver-specific regulation of transcription have been widely reported. as described in methods, the significance of co-occurrence for each a-b tf pair present in the set of sequences was estimated by sampling sequences from the genomic set of sequences until the same number of a and b sites were obtained, and the p-value of frset was computed. in the kidney/liver-specific set we found co-occurrences with p-values <. in one case, the pwm for foxp co-occurred significantly with hnf sites. foxp plays a role in the development of various organs, including liver. in this set of promoters, out of non-overlapping foxp sites co-occurred with hnf sites, yielding a fr value of. given the genomic fr value, we would expect only about foxp sites to co-occur with hnf sites on average. moreover, the distance between foxp sites and hnf sites was biased towards proximal positioning: of the foxp1-hnf site pairs, were separated by less than bps. in addition, visual inspection of the site pairs revealed a preference of the foxp site to be upstream of the hnf site. to illustrate the difference between our approach and approaches based on statistical tests, we calculated co-occurrence p-values using the method of yu et al., and using the method of sudarsanam et al.. the approach by yu et al. evaluates co-occurrences using two p-values, one for co-occurrences, pocc, and one for the bias in distances between pairs of sites, pd. here we focused on pocc, the probability of observing an equal or greater number of co-occurrences, calculated based on the number of sequences in the co-regulated set versus the size of the genome-wide set, the number of co-occurrences between two motifs in the genome-wide set, and the number of co-occurrences in the co-expressed set. the approach by sudarsanam et al. uses a cumulative hypergeometric model to evaluate the significance of the observed number of co-occurrences for a motif pair, by comparing it to the distribution of expected co-occurrences given the number of occurrences of the individual motifs. we applied our fr approach, the pocc approach, and the sudarsanam approach on all sets of co-expressed genes, and compared the results in terms of the over-representation of co-occurring motifs. fig. shows that the distribution of ori p-values for all pwms co-occurring significantly with an over-represented motif is similar to that of all pwms, confirming that the fr approach is not biased by motif over-representation. indeed, the majority of predicted co-occurring motifs are not over-represented. in contrast, the distribution of ori p-values of predicted co-occurring motifs in the top pairs as predicted by pocc, showed a strong bias towards lower ori p-values, indicating that this method is strongly biased by motif over-representation. the fact that with increasing motif over-representation the expected number of co-occurrences modeled by the hypergeometric distribution also increases, makes the approach described by sudarsanam et al. relatively robust against the bias caused by motif over-representation, but less so than the fr measure. however, this method does not use a reference set of sequences during the evaluation of significance, making it the most easily affected of these three approaches by pwm-to-pwm similarities. a relatively high number of co-occurring pairs predicted by the approach by sudarsanam et al., have similar gc content levels, and pairs of motifs with large differences in gc content are relatively rarely predicted to be co-occurring. as an illustration, for the set of promoters of liver- and kidney-specific genes in mouse, the top co-occurrences in terms of pocc were strongly dominated by pwm pairs containing hnf and hnf which were both strongly over-represented in this cluster. in the top motif pairs, involved hnf which was found to have significant pocc values with most other over-represented motifs, such as those for hnf and ikaros. the pair hnf - hnf had the lowest pocc value. however frset was only moderately higher than frgenomic. indeed, only out of hnf sites co-occurred with hnf sites, which were present in out of sequences in this cluster. even though both motifs were over-represented in this cluster, they did not have a strong tendency to be present in the same sequences. the measure described by sudarsanam et al. too, predicted a number of significant co-occurrences involving hnf strikingly, the top motifs predicted to co-occur with hnf motifs have similarly low gc content value as the hnf motif, while this is not the case for the motifs with most significant fr values with regard to hnf. collectively, these results indicate that pocc is more related to co-over-representation than to co-occurrence, and that the measure proposed by sudarsanam et al. is relatively sensitive to gc content similarities. our proposed approach, on the other hand, can find non-over-represented motifs, which are likely to be missed by traditional approaches, and is less influenced by pwm-to-pwm gc content similarity. identification of tfs co-occurring with nf-κb or irf in gene sets having specific expression patterns on tlr stimulation several tfs such as nf-κb, irf, and ap- are known to play roles in gene expression evoked by tlr signaling. however, little is known about involvement of other tfs and how these tfs orchestrate variety of gene expression patterns. to elucidate the diverse gene expression patterns generated by combinations of various tfs, we next applied our method for the analysis on tlr signaling. we used a large scale microarray dataset on gene expression in bone marrow-derived dcs after stimulation with various tlr stimuli. the gene expression values were calculated and hierarchically clustered into clusters. frequency ratio analysis revealed that a wide variety of tfs had significant co-occurrence with each other. firstly, we looked into tfs co-occurring with the binding motifs of nf-κb or irf, since ori analysis showed that these two motifs are significantly over-represented in the clusters, suggesting their key roles in tlr signaling, consistent with former studies. several tfs were found to have significant co-occurrence with nf-κb or irf. nf-κb and irf motifs were found to co-occur significantly with themselves in clusters and clusters and contained genes whose expression levels peaked at about to hours after stimulation, suggesting multiple nf-κb sites in promoters of those genes enable rapid induction. genes in cluster were specific for poly stimulation, and their expression peaked at around hours after stimulation. in addition to various irf family tfs, the pwm m represents stat implying that these genes may be induced secondarily by type i interferon. other well-known motifs found were c-fos, presumably representing the ap- motif. these collectively suggest that our method could identify known tf pairs involved in tlr signaling. we found that in cluster ccaat enhancer binding protein alpha had significant co-occurrence with nf-κb. c/ebp family transcription factors are reported to be involved in tlr signaling-induced gene expression such as cytokine gene expression. importantly, while pocc is relatively low for the pair nf-κb - c/ebpα, up to out of pwms have a pocc < for co-occurrence with the over-represented nf-κb motif in cluster, making pocc not useful for predicting co-occurrences in this case. a gene in the cluster, nos has c/ebp motifs and nf-κb motifs in its promoter, thus the promoter activity is expected to be controlled by c/ebp and/or nf-κb. to check this, we cloned the nos promoter, and its activity upon over-expression of tfs or signaling molecules myd simultaneously with c/ebpα was checked by luciferase assay. the activity of the nos promoter was up-regulated only by over-expression of c/ebpα, rela, or myd, indicating these tfs positively regulate the nos promoter. moreover, when c/ebpα was over-expressed simultaneously with rela or myd luciferase activity increased compared to that on rela or myd over-expression alone. this result indicated that c/ebpα controls the expression of the nos gene, and also suggested that it controls the expression of genes other than nos in clusters synergistic activation of tlr-regulated promoters by nf-κb and c/ebpα the above results prompted us to check if a broader array of promoters is regulated by c/ebpα. we tested whether c/ebpα controls nf-κb-regulated promoters or not. since four pwms in our pwm set represent c/ebp tfs, we picked up promoters having predicted tfbss for nf-κb and one of the four motifs. we also added as a positive control arg which has a predicted nf-κb site and has been reported to be a target of c/ebpβ. five promoters out of tested were up-regulated synergistically by rela and/or myd and c/ebpα, whereas activation of one promoter was suppressed by c/ebpα over-expression. in contrast, the tandem-kb luciferase reporter was suppressed by c/ebpα over-expression, indicating a specific activation of promoters by c/ebpα. these results suggested that c/ebpα synergistically and specifically up-regulates the activity of a set of promoters regulated by nf-κb. taken together, our method successfully identified a pair of transcription factors involved in the immune response.  in this study, we introduced a new measure for regulatory motif co-occurrence, and investigated genome-wide co-occurrence tendencies between pairs of regulatory motifs. our initial results show that some motif pairs have a strong tendency to co-occur, while other pairs have a strong tendency to avoid co-occurrence. however, further investigation showed that these tendencies reflect gc content fluctuations in promoter sequences, rather than a genome-wide level of combinatorial regulation: semi-artificial sequences in which gc content fluctuations were identical to real sequences showed very similar trends. on the other hand, this trend was lost in completely artificial sequences. also, we found that only few pairs of regulatory motifs had exceptionally high or low fr values in the genomic set of sequences as compared to semi-artificial sequences. pairs that did show exceptionally high fr values often involved the tata-box motif, which might reflect the tendency of strictly regulated non-cpg island-associated promoters to contain a tata-box more frequently than other promoter sequences. in addition, we showed that the tendencies were completely different between cpghigh promoters and cpglow promoters. our measure for regulatory motif co-occurrences showed a relatively limited range in cpghigh promoters compared to cpglow promoters. similar observations were made for a set of randomly selected oligonucleotide motifs. these observations reflect a fundamental difference between these two types of promoters. cpg islands have been reported to be associated with ubiquitously expressed genes and housekeeping genes, while genes not associated with cpg islands tend to be tissue-specific or condition-specific genes. on the promoter sequence level too, there are considerable differences: while the promoters of cpg island-associated genes tend to lack typical core promoter elements and tissue-specific tfbss, the promoters of genes not associated with cpg islands tend to contain tata boxes or other core promoter elements and tfbss allowing their precise regulation of expression. the apparent lack of high and low fr values in the genome-wide set of cpghigh promoters might reflect a relatively low need for complex combinatorial regulation, compared to cpglow promoters. on the other hand, for tissue- or condition-specific genes combinatorial regulation might be necessary to ensure spatio-temporal specificity, reflected in the larger range of fr values observed in cpglow promoters. keeping the above observations in mind, in the proposed method, for the set of promoter sequences of interest, the significance of co-occurring pairs was estimated using a random sampling procedure. this approach thus takes into account the genomic tendency of motif pairs of similar structure to appear in the same promoter sequence. furthermore, we considered cpghigh promoters and cpglow promoters as separate cases. a recent study on tfbs analysis has led to a similar recommendation. moreover, we excluded overlapping pairs of sites to avoid bias caused by similarity between motifs. applying our method to a large number of tissue-specific sets of mouse promoters, we could predict a large number of pairs of significantly co-occurring tfbs pairs. one example is the pair hnf - foxp for which we found binding sites to be significantly co-occurring in the promoters of genes with specific expression in liver and kidney. moreover, the tfbss of this pair of tfs showed a tendency to be located proximally to each other, with the foxp tfbss located upstream of the hnf tfbss. importantly, our approach demonstrated improved robustness against biases caused by strongly over-represented motifs in comparison to a previously reported statistics-based method. indeed, the majority of significant interactions we found involved motifs that were not over-represented. this was also the case for c/ebpα binding sites in cluster from the dc expression data. such motifs would thus not be detected by standard over-representation analysis. for one of the significantly co-occurring tf pairs involved in tlr signaling we could verify the predicted combinatorial regulation. we found that c/ebpα co-regulates a set of promoters with nf-κb. co-regulation by nf-κb and c/ebp has previously been reported. lcn and arg were reported as targets of c/ebpβ, confirming that the frequency ratio could predict biologically meaningful tf pairs. moreover, c/ebpα itself has recently been found to control a number of rela-dependent inflammatory promoters, and nf-κb activation synergistically with pu., further supporting our findings on the significance of the nf-κb-c/ebp pair. it would be interesting to check the involvement of other pairs identified as in table in tlr-induced gene expression patterns. there is some room for possible improvements of our approach. first of all, epigenetic factors that might be responsible for tissue- or condition-specific expression should be taken into account. at present, chromatin remodeling data is still limited to a small number of cell types, which makes it difficult to incorporate in our approach. undoubtedly, as the amount of available data increases, there will be a need to incorporate it, resulting in approaches combining both tfbs information and epigenetic information. importantly, while we have limited our analysis to promoter sequences here, the role of distal enhancers in the regulation of transcription is generally accepted. as epigenetic data for various cell types increases, we will become able to apply our method not only to promoter regions but also to enhancers, and investigate potential differences in combinatorial regulation occurring in promoter and enhancers. secondly, since our approach relies on tfbs prediction, which is still known to have a low specificity, further developments in the prediction of tfbss, and additional genome-wide binding data is likely to improve our method as well as other tfbs prediction-based methods. in conclusion, our fr approach circumvents biases which former methodology suffers from, and we could identify some meaningful co-occurring tfbs pairs, one of which was experimentally supported. we believe this approach can help us detect combinatorial interactions between tfs in the regulation of transcription, and we also believe that this sets a basis for future developments in computational identification of combinatorial gene regulation. an online application of our method, which we call regulatory motif combination detector, is available at our website. list of abbreviations used c/ebp: ccaat enhancer binding protein; chip: chromatin immunoprecipitation; dc: dendritic cell; fr: frequency ratio; ori: over-representation index; ppr: pattern-recognition receptor; pwm: position weight matrix; rma: robust multi-array average; tf: transcription factor; tfbs: transcription factor binding site; tlr: toll-like receptor; tss: transcription start site. competing interests the authors declare that they have no competing interests. authors' contributions av and yk conceived of the study. av performed the bioinformatics analysis, and prepared the manuscript. yk carried out experimental validations and prepared the manuscript. dms and sa supervised the project and helped with discussion and interpretation of results and with drafting the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 figure s - workflow of our framework for the detection of co-occurring motifs. the analysis of genome-wide tendencies starts with a set of tfbss, predicted in promoter sequences and a set of pwms. for each pair of motifs, fr values are calculated, and used for further analysis of genome-wide tendencies. the analysis of co-occurrences in sets of co-regulated genes similarly starts with the prediction of tfbss. using these, significantly over-represented tfbss are detected, and for each motif the tendency to co-occur with each of the over-represented motifs is analysed. the significance of the co-occurrences is evaluated using a random sampling approach, sampling sequences from the genomic set of promoters. click here for file additional file 2 supporting text - on the asymmetry of the frequency ration measure. click here for file additional file 3 figure s - genome-wide tendencies of frequency ratios in human promoter sequences. histogram of fr values for all pwm pairs in the genomic set of human promoter sequences. plots of gc content differences as measure of pwm-to-pwm dissimilarity versus fr values, for all promoters, cpghigh promoters, and cpglow promoters. click here for file additional file 4 figure s - tendencies of frequency ratio in semi-artificial and completely artificial sequences. plot of gc content differences as measure of pwm-to-pwm dissimilarity versus fr values in semi-artificial sequences, and completely artificial sequences, semi-artificial cpghigh sequences, and semi-artificial cpglow sequences. click here for file additional file 5 figure s - genome-wide tendencies of frequency ratios for randomly selected 7-mers in human and mouse promoter sequences. plots of gc content differences versus fr values are shown for all human promoters, all mouse promoters, human cpghigh promoters, mouse cpghigh promoters, human cpglow promoters, and mouse cpglow promoters. click here for file additional file 6 figure s - heatmap representation of the average expression values for each of the clusters obtained from the gnf geneatlas mouse data. click here for file additional file 7 table s - summary of main tissues for the clusters obtained from the gnf geneatlas data. click here for file additional file 8 table s - summary of over-represented pwm motifs in tissue-specific sets of mouse promoters click here for file additional file 9 figure s - histogram of the pwm-to-pwm gc content differences of co-occurring motifs predicted by three approaches. co-occurrences predicted by the fr measure are least affected by pwm-to-pwm gc content differences. the distribution of gc content differences of predicted co-occurring pairs of pwms is shown 1) for the pwms found to be significantly co-occurring with an over-represented motif according to fr values, 2) for the pwms found to be co-occurring with an over-represented motif according to pocc, and 3) for the pwms found to be co-occurring with an over-represented motif according to the approach of sudarsanam et al.. for the latter two approaches the pairs with the most significant co-occurrence were used. click here for file additional file 10 figure s - heatmap representation of clusters of tlr-stimulated dc gene expression data referred to in the main text. click here for file additional file 11 table s - summary for the co-occurrences in tissue-specific sets of mouse promoters. click here for file acknowledgements the authors would like to thank all members of the laboratory of systems immunology and the laboratory of host defense for helpful discussions and advice; karlou mar s. amada for help constructing the online tool; dr k. nakai for helpful suggestions; e. kurumatani, a. yoshimura, and y. kimura for technical assistance; m. ogawa, e. kamada and m. kageyama for secretarial assistance. computational time was partly provided by the super computer system at the human genome center, institute of medical sciences, the university of tokyo. funding: this research is supported by the japan society for the promotion of science through the "funding program for world-leading innovative r&d on science and technology," initiated by the council for science and technology policy, and by a kakenhi grant-in-aid for scientific research from the japan society for the promotion of science. this article has been published as part of bmc genomics volume supplement 2012: eleventh international conference on bioinformatics: computational biology. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s 
here we propose a new measure for regulatory motif co-occurrence and a new methodology to systematically identify tf pairs significantly co-occurring in a set of promoter sequences. of these, we experimentally verified that ccaat enhancer binding protein alpha positively regulates its target promoters synergistically with nf-κb. multiple transcription factors are involved in the generation of gene expression patterns, such as tissue-specific gene expression and pleiotropic immune responses. initial analyses suggest that non-cpg promoters have a higher potential for combinatorial regulation than cpg island-associated promoters, and that co-occurrences are strongly influenced by motif similarity. both computational and experimental results indicate that the proposed method can clarify tf interactions that could not be observed by currently available prediction methods. 3- october asia pacific bioinformatics network eleventh international conference on bioinformatics bangkok, thailand binding sites tend to be positioned proximally to each other, suggesting interactions exist between this pair of transcription factors. however, how combinations of tfs orchestrate diverse gene expression patterns is poorly understood. we applied our method to large-scale gene expression data from various tissues, and showed how our measure for motif co-occurrence is not biased by motif over-representation. moreover, the binding sites of several tfs were found to co-occur with nf-κb and irf sites in sets of genes with similar expression patterns in dendritic cells after toll-like receptor stimulation.
 multiple transcription factors are involved in the generation of gene expression patterns, such as tissue-specific gene expression and pleiotropic immune responses. however, how combinations of tfs orchestrate diverse gene expression patterns is poorly understood. here we propose a new measure for regulatory motif co-occurrence and a new methodology to systematically identify tf pairs significantly co-occurring in a set of promoter sequences.  initial analyses suggest that non-cpg promoters have a higher potential for combinatorial regulation than cpg island-associated promoters, and that co-occurrences are strongly influenced by motif similarity. we applied our method to large-scale gene expression data from various tissues, and showed how our measure for motif co-occurrence is not biased by motif over-representation. our method identified, amongst others, the binding motifs of hnf and foxp to be significantly co-occurring in promoters of liver/kidney specific genes. binding sites tend to be positioned proximally to each other, suggesting interactions exist between this pair of transcription factors. moreover, the binding sites of several tfs were found to co-occur with nf-κb and irf sites in sets of genes with similar expression patterns in dendritic cells after toll-like receptor stimulation. of these, we experimentally verified that ccaat enhancer binding protein alpha positively regulates its target promoters synergistically with nf-κb.  both computational and experimental results indicate that the proposed method can clarify tf interactions that could not be observed by currently available prediction methods. 3- october asia pacific bioinformatics network eleventh international conference on bioinformatics bangkok, thailand 