BACKGROUND
the identification of motifs and motif modules is one of the most critical steps to understanding gene regulation.
a motif, often represented by a position weight matrix, is the common pattern of short dna segments bound by a transcription factor .
these dna segments are called transcription factor binding sites .
in high eukaryotes, including humans and mice, it is often the interplay of multiple tfbss from different motifs, instead of from a single motif, that determines the temporal and spatial expression patterns of genes  <cit> .
we define a motif module as a group of several motifs, whose tfbss co-occur in many short dna sequences of one kilobase  long.
we also define cis regulatory modules  as  <dig> kb long sequences containing tfbss of all the motifs of a motif module.
because the binding of tfs to their tfbss plays a pivotal role in controlling gene expression and the dysfunction of these binding sites often results in diseases  <cit> , it is important to identify motifs and motif modules.
many methods are available for the identification of motifs and motif modules.
conventionally, tfbss and motifs were identified on a gene-by-gene basis by experiments such as dnase footprinting  <cit>  and gel-mobility shift assay  <cit> .
it is through these experiments that we understand many basic principles of motifs.
however, such experiments cannot address the challenging problem of identifying motifs and motif modules in the neighborhood of thousands of genes.
thus, many computational methods have been developed  <cit> , which are often based on the following motif properties: overrepresentation, conservation, and clustering.
motif overrepresentation means that tfbss of a motif occur in the non-coding regions of a significant number of genes.
also, tfbss of a motif are often conserved in different species.
finally, motifs are often clustered, with multiple tfbss of different motifs often co-occurring in short dna regions such as crms.
based on these properties, previously developed computational methods have shown some success in identifying motifs and motif modules in a group of putative co-regulated genes as well as on a genome-wide scale.
at the same time, new experimental technologies such as chromatin immunoprecipitation followed by microarray experiments   <cit>  and high-throughput sequencing of immunoprecipitated fragments   <cit>  can provide thousands of short potential tfbs residing regions for computational methods to further identify motifs and motif modules.
although there are many methods for motif and motif module identification, methods that can handle all of the non-coding regions of the human genome are still in great need.
many computational methods were designed to identify novel motifs and motif modules only in the promoter regions of a group of co-regulated genes.
these methods are successful in identifying motifs and motif modules in simple organisms such as yeast.
they are not, however, as successful at identifying motifs and motif modules in higher eukaryotes such as humans.
this is because the tfbss in higher eukaryotes can be several hundred thousand base pairs  upstream, downstream, or in the introns of genes.
to identify tfbss in the long non-coding regions of higher eukaryotes, several methods based on multiple genome alignments were developed  <cit> .
however, current multiple genome alignments may not be able to align tfbss and their orthologous counterparts well.
for instance, multiple genome alignments from several popular methods are significantly different  <cit> .
although chip-chip or chip-seq experiments can narrow down potential tf targeting regions, they are still costly and limited by the availability of high quality antibodies.
thus, novel methods that can systematically identify motif modules in the entire non-coding regions around human genes are needed.
here we describe a computational method that identifies crms and motif modules in the human genome.
unlike the computational methods for promoter regions, our method works on the entire non-coding sequences around human genes.
the non-coding sequences of a human gene include the upstream non-coding sequence until the nearest codon of the 5' adjacent gene, the downstream non-coding sequence until the nearest codon of the 3' adjacent gene, and the intron sequences of the gene itself.
unlike all previous methods, our method measures sequence conservation based on discontiguous sequence similarity  <cit> , which greatly expands the range of the conserved sequences.
our method is also different from the multiple genome alignment based methods, in that we use local alignments, which enables us to identify conserved tfbss and crms that may be "misaligned" in the multiple alignments  <cit> .
by applying the method to all human genes with mouse or rat orthologs in the mouse genome informatics database , we have identified  <dig> motif modules,  <dig> % of which are already supported by various sources of functional evidence.
compared with  <dig> chip-seq experiments, on average, our methods predicted  <dig> % of chip-seq peaks with tfbss of multiple tfs.
our findings also show that tfbss of motifs in many motif modules have preferred distances and orders.
all predicted motif modules are available at http://www.cs.ucf.edu/~xiaoman/module <dig>  we are developing a database, which will enable easier access to these predictions.
RESULTS
a new method for the identification of motif modules and crms
the major obstacle for crm and motif module identification in the entire non-coding regions of the human genome is the large space in which potential tfbss may reside.
for one human gene, the non-coding region can be millions of bp long.
such long non-coding regions for thousands of human genes prevent the direct application of classical statistical methods  <cit>  for crm and motif module identification, due to the convergence time of these methods in the entire human genome.
the large tfbs search space also hinders the direct application of deterministic regular pattern recognition methods because almost every k-mer  is enriched.
note that crms, regulatory units of  <dig> kb long or so that contain multiple tfbss, frequently occur in the human genome and control specific gene expression patterns  <cit> .
therefore, ideally, to understand gene expression patterns, we should identify these regulatory units and consider each unit independently.
that is, we should select short regions that are potential crms from the long non-coding sequences around each gene and then only consider these selected short regions for tfbs identification.
so, how to select potential crms without any idea of the tfbss?
requiring the entire crm regions to be conserved is too restrictive, perhaps only those 6- <dig> bp long tfbss within the crms are conserved.
furthermore, the crm regions may be shuffled in evolution and cannot be aligned with their true orthologous crm regions in multiple genome alignments.
therefore, to identify potential crms, we should develop methods that do not rely on multiple genome alignments.
we developed such a method to systematically identify motif modules and crms in the entire non-coding sequences around human genes .
in brief, we first collected all orthologous gene groups from mgi and defined the non-coding sequences of each gene based on the locations of genes in each species.
each group of orthologous genes contains one human gene and at least one mouse or rat orthologous gene.
the non-coding sequences of a gene include the upstream non-coding sequence until the nearest codon of the 5' adjacent gene, the downstream non-coding sequence until the nearest codon of the 3' adjacent gene, and the intron sequences of the gene itself.
second, from the non-coding regions of every group of orthologous genes, we selected several  <dig> kb long "orthologous" regions based on local alignments and discontiguous sequence similarity .
we call these  <dig> kb long regions blocks.
we used  <dig> kb as the block length cutoff, because  <dig> % of the predicted crms in a genome wide study are shorter than  <dig> kb  <cit> .
third, we defined tfbs candidates by scanning the conserved segments in these blocks using  <dig> vertebrate motifs from the transfac database  <cit> .
finally, we predicted significant motif modules by applying frequent pattern tree  techniques  <cit>  and a poisson clump heuristic  <cit>  to identify frequent instances of motif combinations in these blocks.
in total, we predicted  <dig> motif modules and  <dig> crms.
comparisons with various forms of functional evidence
we compared our predictions with composite element  pairs in the transfac database  <cit>  and pre-defined gene sets in msigdb database  <cit> .
the comparisons show that  <dig> % of the predicted motif modules are supported by functional evidence.
we also compared the predicted tfbss with the readout of  <dig> chip-seq experiments .
with a false discovery rate  < <dig> , the mean of the percentages of the predicted tfbss within the defined chip-seq peaks is  <dig> %, which is comparable with the previous study  <cit> .
we compared the predicted crms with several large-scale crm or enhancer studies as well, which further support our crm and motif module predictions.
in the following, we described each of these comparisons in details to justify our predictions on different levels.
we first compared the predicted motif modules with the ce pairs in transfac  <cit> .
a ce pair contains two different motifs that are experimentally verified to provide combinatorial transcriptional regulation  <cit> .
for the  <dig> vertebrate motifs we used, which can form  <dig> possible motif pairs, there are  <dig> ce pairs in transfac.
for the predicted  <dig> motif modules, there are  <dig> motif pairs, among which  <dig> are ce pairs.
according to the hypergeometric distribution, the p-value of observing at least  <dig> ce motif pairs out of  <dig> predicted motif pairs is  <dig> e- <dig>  several factors, such as the number of orthologous gene groups available, the number of known motifs, and species specific ce pairs, may prevent all  <dig> ce pairs from being included in the predictions.
we next compared the predicted target genes with the collected gene sets in msigdb  <cit> .
the predicted target genes are the closest genes to the predicted crms, which can be in the 5' non-coding regions, the 3' non-coding regions, or the introns of their target genes.
the collected gene sets in msigdb are based on gene annotations, known pathways, gene expression data, mirna target genes and other factors  <cit> .
in table  <dig>  for each category, we listed the number of collected gene sets and the number of motif modules with target genes significantly overlapping with these collected gene sets.
it is surprising to see that the target genes of  <dig>  of the predicted motif modules are also the target genes of mirnas in the pictar gene sets.
mirnas often bind to the 3' of mrnas of their target genes for post-transcriptional control.
the observation that target genes of many motif modules are also target genes of mirnas raises the possibility of coordinated gene regulation by tfs and mirnas at the transcriptional and post-transcriptional level.
with an fdr cutoff of  <dig>  and  <dig> , the target genes of  <dig>  and  <dig>  of the predicted motif modules significantly overlap with at least one gene set in the msigdb database, respectively.
these supported motif modules are most likely biologically meaningful.
in fact, we found literature support for many motif modules.
because of page limit, we only show two examples here.
the first column lists the different sources of functional evidence.
the second column is the number of all gene sets.
the third and fourth columns are the number of motif modules whose target genes significantly overlap with the gene sets.
all of them except kegg pathways are downloaded from msigdb.
for kegg pathways, we directly downloaded from the kegg ftp site.
only the gene sets with at least ten genes are used.
example 1
the motif module "2525420" comprises motifs m <dig> , m <dig> , and m <dig> .
ar is known to interact with oct <dig> physically  <cit> .
oct <dig> is important in lens and olfactory placode development in the mouse  <cit>  and ar is suggested to influence central nervous system development  <cit> .
the function of the tfs in this motif module is consistent with the functions of the target genes.
the  <dig> target genes significantly overlap with the pathway genes in kegg: hsa <dig>  and significantly overlap with the genes annotated with "nervous system development".
axon guidance represents a key stage in the formation of a neuronal network, which supports the functionality of tfs in this motif module.
example 2
the motif module "285236" is composed of motifs m <dig> , m <dig> , m <dig> , and m <dig> .
p <dig> is a well known tumor suppressor protein and involved in various types of cancer.
e2f- <dig> is a member of the e2f family of tfs that plays a crucial role in the control of cell cycle and tumor suppressor protein activities.
ap- <dig> is a tf that regulates proliferation and differentiation in mammalian cells, and is involved in prostate cancer development  <cit> .
krox, also known as egr- <dig>  activates genes that are required for differentiation and mitogenesis, and is over-expressed in prostate cancer  <cit> .
thus, these tfs may work together to play a role in prostate cancer.
consistently, the  <dig> target genes of this motif module significantly overlap with the target genes of mir- <dig>  mir- <dig> induces apoptosis in prostate cancer cells through e2f1-mediated suppression of akt phosphorylation  <cit> , which strongly supports the coordinate gene regulation by this motif module and mir- <dig>
besides the comparisons with the ce pairs in transfac and the gene sets in msigdb, we also compared the predicted tfbss with  <dig> chip-seq experiments .
there are  <dig> motifs corresponding to  <dig> tfs used in these chip-seq experiments.
with an fdr of  <dig> , the mean and the maximal precision of our predictions, defined as the percentage of predicted tfbss within the chip-seq peaks, are  <dig> % and  <dig> %, respectively.
a previous study based on multiple genome alignments had a comparable precision of 3% and 17% in two chip-chip experiments, where the authors selected the predicted crm regions as probes to make the microarray  <cit> .
as in the previous study  <cit> , the precision here is most likely underestimated.
on the other hand, the mean and the maximal recall of our predictions, defined as the percentage of chip-seq peaks with predicted tfbss among all chip-seq peaks with putative binding sites, are  <dig> % and  <dig> %, respectively.
notice that the low recall is mainly due to the fact that chip-seq experiments target tfbs residing regions of individual motifs while our predictions aim to identify motif combinations.
in fact, if we compare the predicted tfbs pairs of two tfs with the chip-seq peaks containing putative tfbss of the two tfs, the mean and the maximal recall is  <dig> % and  <dig> %, respectively, with an fdr of  <dig>  .
for instance, if we consider the individual motifs m <dig>  and m <dig> ,  <dig> % and  <dig> % of the chip-seq peaks are predicted to overlap with crms containing the yy <dig> binding sites and c-myc binding sites, respectively .
if we consider the motif pair ,  <dig> % of the chip-seq peaks with yy <dig> binding sites and c-myc binding sites overlap with the crms containing the two types of binding sites .
we also compared the predicted  <dig> crms with the  <dig> pcrms  <cit> .
these pcrms are computationally predicted crms that are significantly enriched in tfbss for one to five different tfs within  <dig> kb long regions.
to our knowledge, these pcrms are the only available mammalian crms predicted in the entire non-coding sequences around genes based solely on sequence information and known motifs.
identified by using all known vertebrate motifs in the transfac database to scan the genome alignment of human, mouse and rat genomes, these pcrms cover  <dig> human bps .
we found that  <dig> of our predicted  <dig> crms  overlap with the pcrms.
the remaining  <dig> % of crms do not overlap with pcrms, due to the genome rearrangement and/or the fact that genome alignments cannot align tfbss well.
we also found that  <dig> out of  <dig> pcrms  overlap with our predicted crms.
for the remaining  <dig> % of the pcrms,  <dig> % contain tfbss of one or two tfs.
these  <dig> %  of the pcrms may be too weak to be considered as putative functional regions under the schema of the discontiguous sequence similarity.
in addition to the comparison with the pcrms, we also compared our predicted crms with specific enhancers from two recent studies  <cit> .
visel et al.
detected  <dig>   <dig>  and  <dig> p <dig>  peaks in mouse embryonic forebrain, midbrain and limb tissue, respectively  <cit> .
they mapped  <dig> of these sequences to homologous human regions and found that  <dig> of the human regions had enhancer activity.
we found  <dig> out of these  <dig>  enhancers are included in our predicted crms.
narlikar et al.
selected  <dig> human sequences as putative heart enhancers by using known motifs and gene expression data  <cit> .
we found  <dig> out of the  <dig>  predicted enhancers overlap with our predicted crms.
several factors affect the percentage of overlapping with the predicted enhancers in the second study, such as the number of orthologous genes available in our study, the expression data and prior biological knowledge used in the second study, the false positives in each study.
note that published studies on crms often require the occurrence of multiple tfbss in short regions in order to claim these regions to be crms, while we require both the occurrence of multiple tfbss in short regions and the occurrence of the same motif combinations in many short regions.
our requirement of recurrent occurrence of motif combinations may enable us to filter a significant portion of false positives from our predictions.
the order and distance preference of motif modules
the order of the tfbss of different tfs in crms is often important for cooperative tf binding.
by considering the order of a motif pair using the binomial distribution, with an fdr of  <dig> , we found significant orders of  <dig> motif pairs in  <dig> motif modules .
in these motif modules, the median of the percentage of crms with the defined order is  <dig> %.
the order of motif pairs in motif modules supports that these predicted motif modules may be biologically important and the order of tfbss of different tfs may contribute to the efficiency of gene regulation.
these defined orders between motifs will help future motif module identification.
for instance, the motifs m <dig>  and m <dig>  have a preferred order, in which bach <dig> binding sites are on the 5' of the ap- <dig> binding sites.
for each of the  <dig> motif modules containing these two motifs, more than  <dig> % of the crms have this order.
the fact that this order is kept in such a high percentage of the crms of so many motif modules shows that most likely there is an order between the binding sites of the two tfs bach <dig> and ap- <dig>  although there is no literature support this order yet.
such a significant order found in this study may be applied to filter false positive predictions in future crm predictions.
the distance between tfbss within a crm and the distance of the crms to the transcription start  sites often affect the efficiency of gene regulation.
we found that tfbss of  <dig> motif pairs in  <dig> motif modules have specific preferred distance ranges.
for instance, the tfbss of the two motifs m <dig>  and m <dig>  prefer a distance between  <dig> bp and  <dig> bps in our predicted crms.
it has been previously shown that ap2alpha physically interacts with p <dig> to regulate genes in cell growth and metastasis  <cit> .
we also found that the crms of  <dig> motif modules have their distance ranges to the closest transcription start  sites.
for instance, the crms of the motif module "285236" mentioned in example  <dig> have a preferred distance range between  <dig> bp and  <dig> bp to the transcription start sites of their target genes, which further supports the function of this motif module.
furthermore, the crms of  <dig>   <dig>   <dig> motif modules prefer to occur in the upstream, downstream and introns of their target genes, respectively.
the grammar defined by the order and the distance preference of the motif modules will help to discover new crms and filter false positives in future motif identification.
because of the predicted distance preference for a large number of motif modules and the limit of the additional files required by the journal, we can only show some results here.
the detailed information can be found at http://www.cs.ucf.edu/~xiaoman/module <dig>
discussion
we have developed a valuable method to predict motif modules and crms in the entire non-coding regions around human genes.
in total, we have predicted  <dig> crms and  <dig> motif modules,  <dig> % of which are partially supported by various forms of functional evidence with an fdr of  <dig> .
compared with  <dig> chip-seq datasets, on average, with an fdr of  <dig> ,  <dig> % of the predicted tfbss fell into the chip-seq peak regions and  <dig> % of chip-seq peak regions containing multiple tfbss overlap with our predicted crms.
comparing our predictions with the predicted crms in the previous study  <cit> , our method shows at least comparable performance in identifying functional non-coding sequences.
our method is based on local alignments and discontiguous sequence similarity.
unlike coding sequences, tfbss in the non-coding sequences often cannot be aligned well with their corresponding tfbss in other species.
trying to identify crms and motif modules from multiple alignments may miss many crms and motif modules.
moreover, differences among genome alignments generated by different methods affect the results  <cit> .
we find that local alignments are better than multiple alignments for a detailed match of the tfbss and for the identification of functional elements in orthologous sequences.
since often only the segments around tfbss in a crm are conserved, it is natural to consider discontiguous sequence similarity.
to our knowledge, this is the first study that identifies crms in the entire human non-coding regions based on discontiguous sequence similarity.
a significant difference between our predicted crms and those found elsewhere is that our predicted crms share similar motif combinations with many other non-overlapping crm regions.
this recurrence of the same motif combination greatly decreases the chance that the motif instances of the motifs in these combinations co-occur in  <dig> kb long regions by chance.
by comparing the predicted target genes with the gene sets in msigdb, we indeed have shown that  <dig> % of our motif modules are supported by at least one type of evidence.
in future studies, we will include more top  <dig> kb long regions such that the functional elements in the blocks cover a greater percentage of the human genome.
currently, although the total blocks cover about 5% of the human genome, the functional elements within these blocks cover a much smaller percentage of the genome.
we will also incorporate other types of data besides sequence data to identify motif modules and crms.
for instance, gene expression data and gene tissue specificity will definitely help to filter false positive predictions.
moreover, incorporating information from recent high throughput enhancer studies  <cit>  should further improve the developed method.
with these additions, we expect to identify more significant motif modules.
moreover, the sensitivity of all the predicted motif modules could be further improved.
CONCLUSIONS
our work provides a large-scale prediction of motif modules in mammals, and will facilitate the understanding of gene regulation in a systematic way.

 the identification of motifs and motif modules is one of the most critical steps to understanding gene regulation. a motif, often represented by a position weight matrix, is the common pattern of short dna segments bound by a transcription factor. these dna segments are called transcription factor binding sites. in high eukaryotes, including humans and mice, it is often the interplay of multiple tfbss from different motifs, instead of from a single motif, that determines the temporal and spatial expression patterns of genes. we define a motif module as a group of several motifs, whose tfbss co-occur in many short dna sequences of one kilobase long. we also define cis regulatory modules as kb long sequences containing tfbss of all the motifs of a motif module. because the binding of tfs to their tfbss plays a pivotal role in controlling gene expression and the dysfunction of these binding sites often results in diseases, it is important to identify motifs and motif modules. many methods are available for the identification of motifs and motif modules. conventionally, tfbss and motifs were identified on a gene-by-gene basis by experiments such as dnase footprinting and gel-mobility shift assay. it is through these experiments that we understand many basic principles of motifs. however, such experiments cannot address the challenging problem of identifying motifs and motif modules in the neighborhood of thousands of genes. thus, many computational methods have been developed, which are often based on the following motif properties: overrepresentation, conservation, and clustering. motif overrepresentation means that tfbss of a motif occur in the non-coding regions of a significant number of genes. also, tfbss of a motif are often conserved in different species. finally, motifs are often clustered, with multiple tfbss of different motifs often co-occurring in short dna regions such as crms. based on these properties, previously developed computational methods have shown some success in identifying motifs and motif modules in a group of putative co-regulated genes as well as on a genome-wide scale. at the same time, new experimental technologies such as chromatin immunoprecipitation followed by microarray experiments and high-throughput sequencing of immunoprecipitated fragments can provide thousands of short potential tfbs residing regions for computational methods to further identify motifs and motif modules. although there are many methods for motif and motif module identification, methods that can handle all of the non-coding regions of the human genome are still in great need. many computational methods were designed to identify novel motifs and motif modules only in the promoter regions of a group of co-regulated genes. these methods are successful in identifying motifs and motif modules in simple organisms such as yeast. they are not, however, as successful at identifying motifs and motif modules in higher eukaryotes such as humans. this is because the tfbss in higher eukaryotes can be several hundred thousand base pairs upstream, downstream, or in the introns of genes. to identify tfbss in the long non-coding regions of higher eukaryotes, several methods based on multiple genome alignments were developed. however, current multiple genome alignments may not be able to align tfbss and their orthologous counterparts well. for instance, multiple genome alignments from several popular methods are significantly different. although chip-chip or chip-seq experiments can narrow down potential tf targeting regions, they are still costly and limited by the availability of high quality antibodies. thus, novel methods that can systematically identify motif modules in the entire non-coding regions around human genes are needed. here we describe a computational method that identifies crms and motif modules in the human genome. unlike the computational methods for promoter regions, our method works on the entire non-coding sequences around human genes. the non-coding sequences of a human gene include the upstream non-coding sequence until the nearest codon of the 5' adjacent gene, the downstream non-coding sequence until the nearest codon of the 3' adjacent gene, and the intron sequences of the gene itself. unlike all previous methods, our method measures sequence conservation based on discontiguous sequence similarity, which greatly expands the range of the conserved sequences. our method is also different from the multiple genome alignment based methods, in that we use local alignments, which enables us to identify conserved tfbss and crms that may be "misaligned" in the multiple alignments. by applying the method to all human genes with mouse or rat orthologs in the mouse genome informatics database, we have identified motif modules, % of which are already supported by various sources of functional evidence. compared with chip-seq experiments, on average, our methods predicted % of chip-seq peaks with tfbss of multiple tfs. our findings also show that tfbss of motifs in many motif modules have preferred distances and orders. all predicted motif modules are available at http://www.cs.ucf.edu/~xiaoman/module we are developing a database, which will enable easier access to these predictions.  a new method for the identification of motif modules and crms the major obstacle for crm and motif module identification in the entire non-coding regions of the human genome is the large space in which potential tfbss may reside. for one human gene, the non-coding region can be millions of bp long. such long non-coding regions for thousands of human genes prevent the direct application of classical statistical methods for crm and motif module identification, due to the convergence time of these methods in the entire human genome. the large tfbs search space also hinders the direct application of deterministic regular pattern recognition methods because almost every k-mer is enriched. note that crms, regulatory units of kb long or so that contain multiple tfbss, frequently occur in the human genome and control specific gene expression patterns. therefore, ideally, to understand gene expression patterns, we should identify these regulatory units and consider each unit independently. that is, we should select short regions that are potential crms from the long non-coding sequences around each gene and then only consider these selected short regions for tfbs identification. so, how to select potential crms without any idea of the tfbss? requiring the entire crm regions to be conserved is too restrictive, perhaps only those 6- bp long tfbss within the crms are conserved. furthermore, the crm regions may be shuffled in evolution and cannot be aligned with their true orthologous crm regions in multiple genome alignments. therefore, to identify potential crms, we should develop methods that do not rely on multiple genome alignments. we developed such a method to systematically identify motif modules and crms in the entire non-coding sequences around human genes. in brief, we first collected all orthologous gene groups from mgi and defined the non-coding sequences of each gene based on the locations of genes in each species. each group of orthologous genes contains one human gene and at least one mouse or rat orthologous gene. the non-coding sequences of a gene include the upstream non-coding sequence until the nearest codon of the 5' adjacent gene, the downstream non-coding sequence until the nearest codon of the 3' adjacent gene, and the intron sequences of the gene itself. second, from the non-coding regions of every group of orthologous genes, we selected several kb long "orthologous" regions based on local alignments and discontiguous sequence similarity. we call these kb long regions blocks. we used kb as the block length cutoff, because % of the predicted crms in a genome wide study are shorter than kb. third, we defined tfbs candidates by scanning the conserved segments in these blocks using vertebrate motifs from the transfac database. finally, we predicted significant motif modules by applying frequent pattern tree techniques and a poisson clump heuristic to identify frequent instances of motif combinations in these blocks. in total, we predicted motif modules and crms. comparisons with various forms of functional evidence we compared our predictions with composite element pairs in the transfac database and pre-defined gene sets in msigdb database. the comparisons show that % of the predicted motif modules are supported by functional evidence. we also compared the predicted tfbss with the readout of chip-seq experiments. with a false discovery rate <, the mean of the percentages of the predicted tfbss within the defined chip-seq peaks is %, which is comparable with the previous study. we compared the predicted crms with several large-scale crm or enhancer studies as well, which further support our crm and motif module predictions. in the following, we described each of these comparisons in details to justify our predictions on different levels. we first compared the predicted motif modules with the ce pairs in transfac. a ce pair contains two different motifs that are experimentally verified to provide combinatorial transcriptional regulation. for the vertebrate motifs we used, which can form possible motif pairs, there are ce pairs in transfac. for the predicted motif modules, there are motif pairs, among which are ce pairs. according to the hypergeometric distribution, the p-value of observing at least ce motif pairs out of predicted motif pairs is e- several factors, such as the number of orthologous gene groups available, the number of known motifs, and species specific ce pairs, may prevent all ce pairs from being included in the predictions. we next compared the predicted target genes with the collected gene sets in msigdb. the predicted target genes are the closest genes to the predicted crms, which can be in the 5' non-coding regions, the 3' non-coding regions, or the introns of their target genes. the collected gene sets in msigdb are based on gene annotations, known pathways, gene expression data, mirna target genes and other factors. in table for each category, we listed the number of collected gene sets and the number of motif modules with target genes significantly overlapping with these collected gene sets. it is surprising to see that the target genes of of the predicted motif modules are also the target genes of mirnas in the pictar gene sets. mirnas often bind to the 3' of mrnas of their target genes for post-transcriptional control. the observation that target genes of many motif modules are also target genes of mirnas raises the possibility of coordinated gene regulation by tfs and mirnas at the transcriptional and post-transcriptional level. with an fdr cutoff of and, the target genes of and of the predicted motif modules significantly overlap with at least one gene set in the msigdb database, respectively. these supported motif modules are most likely biologically meaningful. in fact, we found literature support for many motif modules. because of page limit, we only show two examples here. the first column lists the different sources of functional evidence. the second column is the number of all gene sets. the third and fourth columns are the number of motif modules whose target genes significantly overlap with the gene sets. all of them except kegg pathways are downloaded from msigdb. for kegg pathways, we directly downloaded from the kegg ftp site. only the gene sets with at least ten genes are used. example 1 the motif module "2525420" comprises motifs m, m, and m. ar is known to interact with oct physically. oct is important in lens and olfactory placode development in the mouse and ar is suggested to influence central nervous system development. the function of the tfs in this motif module is consistent with the functions of the target genes. the target genes significantly overlap with the pathway genes in kegg: hsa and significantly overlap with the genes annotated with "nervous system development". axon guidance represents a key stage in the formation of a neuronal network, which supports the functionality of tfs in this motif module. example 2 the motif module "285236" is composed of motifs m, m, m, and m. p is a well known tumor suppressor protein and involved in various types of cancer. e2f- is a member of the e2f family of tfs that plays a crucial role in the control of cell cycle and tumor suppressor protein activities. ap- is a tf that regulates proliferation and differentiation in mammalian cells, and is involved in prostate cancer development. krox, also known as egr- activates genes that are required for differentiation and mitogenesis, and is over-expressed in prostate cancer. thus, these tfs may work together to play a role in prostate cancer. consistently, the target genes of this motif module significantly overlap with the target genes of mir- mir- induces apoptosis in prostate cancer cells through e2f1-mediated suppression of akt phosphorylation, which strongly supports the coordinate gene regulation by this motif module and mir- besides the comparisons with the ce pairs in transfac and the gene sets in msigdb, we also compared the predicted tfbss with chip-seq experiments. there are motifs corresponding to tfs used in these chip-seq experiments. with an fdr of, the mean and the maximal precision of our predictions, defined as the percentage of predicted tfbss within the chip-seq peaks, are % and %, respectively. a previous study based on multiple genome alignments had a comparable precision of 3% and 17% in two chip-chip experiments, where the authors selected the predicted crm regions as probes to make the microarray. as in the previous study, the precision here is most likely underestimated. on the other hand, the mean and the maximal recall of our predictions, defined as the percentage of chip-seq peaks with predicted tfbss among all chip-seq peaks with putative binding sites, are % and %, respectively. notice that the low recall is mainly due to the fact that chip-seq experiments target tfbs residing regions of individual motifs while our predictions aim to identify motif combinations. in fact, if we compare the predicted tfbs pairs of two tfs with the chip-seq peaks containing putative tfbss of the two tfs, the mean and the maximal recall is % and %, respectively, with an fdr of. for instance, if we consider the individual motifs m and m, % and % of the chip-seq peaks are predicted to overlap with crms containing the yy binding sites and c-myc binding sites, respectively. if we consider the motif pair, % of the chip-seq peaks with yy binding sites and c-myc binding sites overlap with the crms containing the two types of binding sites. we also compared the predicted crms with the pcrms. these pcrms are computationally predicted crms that are significantly enriched in tfbss for one to five different tfs within kb long regions. to our knowledge, these pcrms are the only available mammalian crms predicted in the entire non-coding sequences around genes based solely on sequence information and known motifs. identified by using all known vertebrate motifs in the transfac database to scan the genome alignment of human, mouse and rat genomes, these pcrms cover human bps. we found that of our predicted crms overlap with the pcrms. the remaining % of crms do not overlap with pcrms, due to the genome rearrangement and/or the fact that genome alignments cannot align tfbss well. we also found that out of pcrms overlap with our predicted crms. for the remaining % of the pcrms, % contain tfbss of one or two tfs. these % of the pcrms may be too weak to be considered as putative functional regions under the schema of the discontiguous sequence similarity. in addition to the comparison with the pcrms, we also compared our predicted crms with specific enhancers from two recent studies. visel et al. detected and p peaks in mouse embryonic forebrain, midbrain and limb tissue, respectively. they mapped of these sequences to homologous human regions and found that of the human regions had enhancer activity. we found out of these enhancers are included in our predicted crms. narlikar et al. selected human sequences as putative heart enhancers by using known motifs and gene expression data. we found out of the predicted enhancers overlap with our predicted crms. several factors affect the percentage of overlapping with the predicted enhancers in the second study, such as the number of orthologous genes available in our study, the expression data and prior biological knowledge used in the second study, the false positives in each study. note that published studies on crms often require the occurrence of multiple tfbss in short regions in order to claim these regions to be crms, while we require both the occurrence of multiple tfbss in short regions and the occurrence of the same motif combinations in many short regions. our requirement of recurrent occurrence of motif combinations may enable us to filter a significant portion of false positives from our predictions. the order and distance preference of motif modules the order of the tfbss of different tfs in crms is often important for cooperative tf binding. by considering the order of a motif pair using the binomial distribution, with an fdr of, we found significant orders of motif pairs in motif modules. in these motif modules, the median of the percentage of crms with the defined order is %. the order of motif pairs in motif modules supports that these predicted motif modules may be biologically important and the order of tfbss of different tfs may contribute to the efficiency of gene regulation. these defined orders between motifs will help future motif module identification. for instance, the motifs m and m have a preferred order, in which bach binding sites are on the 5' of the ap- binding sites. for each of the motif modules containing these two motifs, more than % of the crms have this order. the fact that this order is kept in such a high percentage of the crms of so many motif modules shows that most likely there is an order between the binding sites of the two tfs bach and ap- although there is no literature support this order yet. such a significant order found in this study may be applied to filter false positive predictions in future crm predictions. the distance between tfbss within a crm and the distance of the crms to the transcription start sites often affect the efficiency of gene regulation. we found that tfbss of motif pairs in motif modules have specific preferred distance ranges. for instance, the tfbss of the two motifs m and m prefer a distance between bp and bps in our predicted crms. it has been previously shown that ap2alpha physically interacts with p to regulate genes in cell growth and metastasis. we also found that the crms of motif modules have their distance ranges to the closest transcription start sites. for instance, the crms of the motif module "285236" mentioned in example have a preferred distance range between bp and bp to the transcription start sites of their target genes, which further supports the function of this motif module. furthermore, the crms of motif modules prefer to occur in the upstream, downstream and introns of their target genes, respectively. the grammar defined by the order and the distance preference of the motif modules will help to discover new crms and filter false positives in future motif identification. because of the predicted distance preference for a large number of motif modules and the limit of the additional files required by the journal, we can only show some results here. the detailed information can be found at http://www.cs.ucf.edu/~xiaoman/module discussion we have developed a valuable method to predict motif modules and crms in the entire non-coding regions around human genes. in total, we have predicted crms and motif modules, % of which are partially supported by various forms of functional evidence with an fdr of. compared with chip-seq datasets, on average, with an fdr of, % of the predicted tfbss fell into the chip-seq peak regions and % of chip-seq peak regions containing multiple tfbss overlap with our predicted crms. comparing our predictions with the predicted crms in the previous study, our method shows at least comparable performance in identifying functional non-coding sequences. our method is based on local alignments and discontiguous sequence similarity. unlike coding sequences, tfbss in the non-coding sequences often cannot be aligned well with their corresponding tfbss in other species. trying to identify crms and motif modules from multiple alignments may miss many crms and motif modules. moreover, differences among genome alignments generated by different methods affect the results. we find that local alignments are better than multiple alignments for a detailed match of the tfbss and for the identification of functional elements in orthologous sequences. since often only the segments around tfbss in a crm are conserved, it is natural to consider discontiguous sequence similarity. to our knowledge, this is the first study that identifies crms in the entire human non-coding regions based on discontiguous sequence similarity. a significant difference between our predicted crms and those found elsewhere is that our predicted crms share similar motif combinations with many other non-overlapping crm regions. this recurrence of the same motif combination greatly decreases the chance that the motif instances of the motifs in these combinations co-occur in kb long regions by chance. by comparing the predicted target genes with the gene sets in msigdb, we indeed have shown that % of our motif modules are supported by at least one type of evidence. in future studies, we will include more top kb long regions such that the functional elements in the blocks cover a greater percentage of the human genome. currently, although the total blocks cover about 5% of the human genome, the functional elements within these blocks cover a much smaller percentage of the genome. we will also incorporate other types of data besides sequence data to identify motif modules and crms. for instance, gene expression data and gene tissue specificity will definitely help to filter false positive predictions. moreover, incorporating information from recent high throughput enhancer studies should further improve the developed method. with these additions, we expect to identify more significant motif modules. moreover, the sensitivity of all the predicted motif modules could be further improved.  our work provides a large-scale prediction of motif modules in mammals, and will facilitate the understanding of gene regulation in a systematic way. 
in total, with a false discovery rate cutoff of, we predicted, motif modules, % of which are supported by various forms of functional evidence. to identify motif modules systematically, we developed a computational method for the entire non-coding regions around human genes that does not rely upon the use of multiple genome alignments. finally, a frequent pattern mining technique was applied to identify motif modules within these blocks. next, we scanned the conserved segments in these blocks using known motifs in the transfac database. compared with experimental data from chip-seq experiments, on average, our methods predicted % of the chip-seq peaks with tfbss of multiple tfs. our work provides a large-scale prediction of motif modules in mammals, which will facilitate the understanding of gene regulation in a systematic way. however, the promoter regions only account for a limited number of locations where transcription factor binding sites can occur, and multiple genome alignments often cannot align binding sites with their true counterparts because of the short and degenerative nature of these transcription factor binding sites. our findings also show that many motif modules have distance preference and order preference among the motifs, which further supports the functionality of these predictions. current approaches to identifying motif modules are often restricted to searches within promoter regions or rely on multiple genome alignments. first, we selected orthologous dna blocks approximately 1-kilobase in length based on discontiguous sequence similarity. the identification of motif modules, groups of multiple motifs frequently occurring in dna sequences, is one of the most important tasks necessary for annotating the human genome.
 the identification of motif modules, groups of multiple motifs frequently occurring in dna sequences, is one of the most important tasks necessary for annotating the human genome. current approaches to identifying motif modules are often restricted to searches within promoter regions or rely on multiple genome alignments. however, the promoter regions only account for a limited number of locations where transcription factor binding sites can occur, and multiple genome alignments often cannot align binding sites with their true counterparts because of the short and degenerative nature of these transcription factor binding sites.  to identify motif modules systematically, we developed a computational method for the entire non-coding regions around human genes that does not rely upon the use of multiple genome alignments. first, we selected orthologous dna blocks approximately 1-kilobase in length based on discontiguous sequence similarity. next, we scanned the conserved segments in these blocks using known motifs in the transfac database. finally, a frequent pattern mining technique was applied to identify motif modules within these blocks. in total, with a false discovery rate cutoff of, we predicted, motif modules, % of which are supported by various forms of functional evidence. compared with experimental data from chip-seq experiments, on average, our methods predicted % of the chip-seq peaks with tfbss of multiple tfs. our findings also show that many motif modules have distance preference and order preference among the motifs, which further supports the functionality of these predictions.  our work provides a large-scale prediction of motif modules in mammals, which will facilitate the understanding of gene regulation in a systematic way. 