BACKGROUND
to understand how gene expression is coordinated to produce hundreds of cell phenotypes from an identical complement of genes is a principal challenge in mammalian genome research.
a commonly suggested model for terminal differentiation in metazoans is that the core features of the cellular phenotype are mediated by a set of genes that is regulated as a gene battery, i.e.
a set of functionally coupled genes that are activated by similar cis- and trans-acting regulators  <cit> .
although the gene battery is an idealized concept, concrete examples of gene battery-like regulation have been found in for example muscle subtypes  <cit> , megakaryocytes  <cit> , the epidermis  <cit>  and lymphocytes  <cit> .
a key step in the elucidation of gene battery-like regulation is to detect and functionally test the cis regulatory elements that mediate the co-regulation.
a number of computation-based methods have been proposed to do this.
in micro-organisms, computational methods have proven useful to detect modules of co-regulated genes  <cit> .
in mammals, predictive models based on assumed co-regulation at the cis level have been constructed for liver- and skeletal muscle-selective gene regulation  <cit> , and general tools have been developed for the regulatory analysis of co-expressed genes  <cit> .
the aim of this work is to screen the mouse and human genomes and transcriptomes for instances where sharing of cis-regulatory sequences is statistically coupled to conserved co-expression of genes, i.e.
cases that fall within or near the idealized gene battery concept.
another aim is to critically investigate technical parameters in order to maximize the sensitivity by which co-regulation of co-expressed genes can be detected.
in a screen that covered ~ <dig> per cent of all protein-coding genes according to the latest ensembl annotation, we identified  <dig> co-expressed gene clusters with  <dig> cases of statistically supported sharing of cis-regulatory sequence motifs.
the predictive value of the assignment of transcription factor binding sites was experimentally evaluated on ebf binding sites in a set of b-cell expressed gene clusters.
the predicted cases of co-regulation include several previously known prototype examples of tissue specific regulation, but also novel predictions.
all data are made available to the research community in the form of an internet resource that may serve as a starting point for further analysis.
RESULTS
the analysis was based on the assumption that homologous genes in mouse and human are equivalent in most aspects of regulation and function.
in particular, we assumed that the transcriptional regulation is conserved for orthologous genes.
for example, the mouse gene myh <dig> and the human gene myh <dig> are assumed to share expression pattern and to share important cis-regulatory sequences.
below, the term 'ortholog pair' will be applied as a two-species equivalent of 'gene', for which expression and sequence data where retrieved for both mouse and human.
ensembl gene annotations  <cit>  were applied thoughout the analysis.
co-expressed gene sets were defined from a compendium of mouse and human expression data and tested for functional coupling
previous results by our group and others have shown that statistical analysis of gene expression profiles in a large compendium of expression data can predict targets of differentiation processes, and identify functionally coupled genes  <cit> .
in the following analysis, we specifically focused on a compendium derived from the recently completed novartis expression atlas   <cit> .
these data contain transcription profiles for  <dig> mouse and human tissues generated by hybridization on customized affymetrix chips, and cover a large fraction of the mouse and human protein-coding genes.
sequence annotation of the novartis probes linked the mouse and human data to approximately  <dig>  unique ensembl gene identities in each species .
between the two datasets,  <dig> non-redundant ortholog pairs could be identified by linkage of reciprocal ensembl homology assignments .
in later steps, we excluded genes for which regulatory sequence could not be extracted , leaving  <dig> ortholog pairs for clustering .
the final dataset included ~40% of the mouse and human ensembl annotated genes.
numbers in the 'mouse' and 'human' columns signify the number of unique ensembl identifiers in each respective species.
numbers in column 'both' signify ortholog pairs of ensembl entries.
the overlap between the nonredundant sequence database  and the nonredundant expression dataset  was  <dig> ortholog pairs.
clustering
we clustered the mouse/human ortholog pairs based on their expression profiles across the  <dig> mouse and human tissues .
we computed clusters at cut-off levels ranging from pearson's correlation coefficient   <dig>  to  <dig> .
at the lowest applied cut-off, 57% of all ortholog pairs in the data were members of a cluster.
the cluster sizes were distributed in a skewed manner, with a predominant formation of small clusters .
all analyses hereafter were performed at a pcc =  <dig>  cut-off, which generated  <dig> clusters with  <dig> ortholog pairs.
this relatively stringent cut-off was chosen to reduce the number of non-relevant genes in the clusters.
a higher cut-off did not seem reasonable given the noise-level of the microarray experiments .
assessment of functional linkage
according to the definition, a gene battery should encode functionally linked proteins.
we used gene ontology  terms, protein-protein interactions , and manual curation to assess functional linkage within our expression clusters.
first, we investigated the relationship between expression profile similarity of two ortholog pairs and their relative probability to share a functional annotation term or to encode interacting proteins .
there was a consistent correlation between co-expression and the relative probability for two genes to share a go-term or to encode interacting proteins .
second, we studied the statistical over-representation of go terms and interacting proteins inside the clusters .
30/ <dig> clusters with ten or more ortholog pairs contained at least one over-represented go term.
the proportion of small clusters with over-represented go terms was lower, which reflects a lack of statistical power in small clusters.
genes encoding interacting protein pairs were also over-represented inside clusters.
<dig> cases of protein-protein interaction between two genes in the same cluster were found.
in contrast,  <dig> simulations on permuted data revealed a median of  <dig> interactions .
the bind data contained only  <dig> interactions that could be mapped to the dataset, which explains the seemingly low number of  <dig> interactions.
finally, the clusters were annotated by manual curation.
six examples are shown in figures  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>  and clusters with predicted regulators are listed in table  <dig>  for a full overview, see the web supplement .
several clusters clearly represent gene sets that mediate specialized features of different cell types, including smooth muscle specific genes , b lymphocyte-specific genes , and genes selectively expressed in the testis .
further, we detected clusters that were related to cellular processes or organelles, including endoplasmatic reticulum , protein synthesis , and modulators of transcriptional regulation .
a grand majority of the clusters were defined by peaks at different amplitudes in several tissues.
as an example, the endoplasmatic reticulum cluster  was defined by a highly variant profile with strong expression in, for example, exocrine glands.
generally, the cluster profiles were conserved between species, in the sense that clusters were defined by peaks in the same organs.
this effect was more pronounced for clusters with expression in a single organ, such as the testis .
over-represented motifs arranged by cluster number.
fdr column: false discovery rate .
motifs are are shown both by their numerical identifiers  and by their annotation .
in cases where a pfm is a composite based on more than one source, the components are given separated by commas.
the data where generated from the pcc =  <dig>  clustering,  <dig> kb sequence database, at 90% phylogenetic conservation.
in combination, the go term enrichment, the protein-protein interaction, and the manual curation convincingly show that clustered genes are functionally linked.
regulatory dna and descriptions of transcription factor binding sites were extracted and pre-processed
in the next part of the analysis, individual ortholog pairs were scored for transcription factor binding sites.
binding motifs were represented in the form of position frequency matrices .
based on a fixed amount of upstream dna sequence in each ortholog pair, a statistical score was computed to predict the potential for a site in the sequence to bind the factor corresponding to a pfm .
we extracted upstream + intronic sequence from the ensembl database, in amounts of  <dig>   <dig> or  <dig> kb per gene .
a filter was applied that removed ortholog pairs for which the transcription start differed between the two species .
filtering was successful for  <dig> unique ortholog pairs .
in a subsequent filter, dna sequence that was not conserved between mouse and human was removed, so-called phylogenetic footprinting.
phylogenetic footprinting was applied at different stringency, to allow the following optimization of the protocol .
furthermore, all exon sequence was removed from the analysis .
finally, the sequences were matched with the expression data based on annotation, the overlap being  <dig> ortholog pairs .
<dig> vertebrate pfms were downloaded from the transfac and jaspar databases  <cit> .
since the databases appeared to contain redundant or equivalent entries, highly similar pfms were grouped and merged using single linkage hierarchical clustering and a pfm distance measure defined in  <cit> , which reduced the number of pfms from  <dig> to  <dig> .
this step reduced the redundancy, but did not merge all identically annotated pfms  pfms in table  <dig>  cluster 40).
design of a predictive scoring system for transcription factor binding
after the retrieval and preprocessing of both sequences and pfms, all individual sequences were tested for pfm matches using the mast software  <cit> , a software for identifying single or multiple motifs in sequences.
mast was set to compute one p-value for each pfm with respect to each sequence .
based on the p-values obtained from the mast software, a composite score was defined as the product of the p-value in the mouse and human sequences of an ortholog pair .
a composite score close to  <dig> indicates that both the mouse and human promoter sequence in the ortholog pair contains sequence elements that are in very good agreement with a certain pfm.
to address the biological validity of the mast composite scores within the context of a set of co-expressed genes, we screened  <dig> ortholog pairs present in b-cell expressed clusters for individual ebf sites.
in all,  <dig> individual ebf sites in  <dig> different ortholog pairs were detected .
to test the functionality  of these sites, we examined the potential of  <dig> basepair duplex oligo-nucleotides spanning the presumed sites to compete for protein binding in electrophoretic mobility shift assays .
ebf binding capacity was assayed using nuclear extracts from the pre-b cell line 40-ei, a labelled mouse mb- <dig> promoter high affinity ebf site  <cit> , and competitor oligo-nucleotides covering the new potential binding sites.
in the absence of competitor oligo-nucleotide, a prominent dna/protein complex  could be detected whereas this complex was undetectable after the addition of a  <dig> fold molar excess of the unlabeled binding site .
the identity of the protein and the specificity of the binding were verified by competition with a point mutated mb-1/ebf site and by the inclusion of an ebf specific antibody into the reaction mixture .
the point mutated ebf binding site was unable to abolish complex formation even in a 1000-fold molar excess , indicating that we detect specific protein dna interactions with this experimental set up.
<dig> out of the  <dig> new binding sites competed for complex formation when added in a 300- or 1000-fold molar excess, and thus have the ability to bind ebf in vitro .
we conclude that the large majority of binding sites were able to functionally interact with the predicted protein, and that the composite score in principle detected factor binding.
it should however be emphasized that the quality of predictions is dependent on the quality of the binding site descriptions, and the result does not necessarily imply that other predicted factors bind.
a statistical procedure was used to detect enriched motifs in the clusters
to test whether the identified clusters represent potential gene batteries, i.e.
contain shared cis-regulatory elements, we designed a procedure to detect significant over-representation of orthologs that match a pfm inside a cluster.
the procedure is based on a modification of fisher's exact test, which tests for dependency between two events   <cit> .
we introduced a procedure to optimize the composite score thresholds for individual pfms.
in brief, we selected the threshold that gave the lowest fisher test p-value in any one cluster.
this was based on the assumption that non-random distribution of detections over clusters reflects biological function, as has been proposed in  <cit> .
the fisher test p-values with the optimized thresholds are hereafter termed p-scores.
the tests of multiple detection thresholds, multiple clusters and  <dig> pfms led to a need to compensate for mass testing.
this was done by estimating false discovery rates  based on simulations on randomized data.
in brief, we compared the outcome when using permuted and when using observed data at different p-scores, and defined the fdr as the ratio between the two .
this procedure allowed us to choose a significance threshold with a controlled expected number of spurious detections.
simulations were repeated  <dig> times.
since optimization of detection thresholds were repeated in each simulation round, no bias in disfavour of the control case was introduced.
the detection of over-represented pfms was affected by dna amount and masking, but not affected by gene paralogy
using the described algorithm, we analyzed the number of times a motif was over-represented when using different amounts of dna sequence per gene  and different stringency in the phylogenetic footprinting , >60%, >70%, >80% or >90% identity).
the analysis was performed at 10% false discovery rate and with the clustering obtained at pcc =  <dig>  clustering cut-off.
at all dna amounts, higher stringency phylogenetic footprinting appeared to be beneficial, and between dna amounts,  <dig> kb and  <dig> kb compared favourably over  <dig> kb .
we conclude that optimal results are obtained when using a limited amount of sequence per gene .
a potential confounding factor in the analysis is that similarity in upstream sequence may be attributable to factors other than shared cis-regulatory elements.
the most important such factor is likely to be gene paralogy, since over-represented motifs might simply represent matches to non-functional  sequence in a co-expressed gene family.
in total,  <dig> cases of pairwise paralogy were detected in the dataset of  <dig> ortholog pairs using blast.
all genes with transcripts matches at a blast e-value less than  <dig>  were defined as potential paralogs.
of the  <dig> ortholog pairs that clustered,  <dig> genes had a paralog inside the same cluster.
this was clearly an over-representation since filling clusters with random genes from the dataset produced an average of  <dig> genes with a paralog inside a cluster .
in a side-by-side comparison, we therefore analyzed the effect of removing paralogs from a cluster before testing for over-represented motifs , as opposed to not removing paralogs but the same number of randomly selected genes.
the results showed no increase in the number of over-represented motifs when allowing paralogs in the same cluster .
<dig> over-represented motifs were observed in  <dig> clusters
the final analysis was performed with  <dig> kb sequence and 80% phylogenetic footprinting .
again, the clustering generated at pcc =  <dig>  was chosen for analysis.
the over-representation algorithm was run in the above cases, using  <dig> simulations to estimate false discovery rate thresholds.
over-represented motifs with fdr:s less than  <dig> % and 10% were recorded.
key features of these results are presented in tables  <dig> and  <dig>  and the complete results are available by web browser .
criteria for inclusion in the different groups: i. existing evidence of co-regulation of cluster member by the factor, and/or a mutant mouse phenotype that affects differentiation of the cell type.
ii.
evidence of the factor acting on a limited number of genes in the tissue.
evidence of the transcription factor being expressed in the tissue.
iii.
no or limited evidence in the literature.
* = detected at 20% fdr.
a predictive screen performed with these settings associated  <dig> motifs to co-expression in a total of  <dig> clusters at fdr thresholds  <dig> % and 10% .
there was a clear tendency for motifs to be detected in most large clusters, and for smaller clusters to lack over-represented motifs.
this can be directly explained by the fact that even a high degree of motif presence in a small cluster can be attributed to spurious detections, and that the fisher test accounts for that.
the over-represented motifs were validated against the literature, and fell into three main categories: 1) previously described cases of gene-battery-like regulation, 2)previously unreported cases of gene battery-like regulation with some support in a limited number of genes, and 3) hypothetical cases of gene battery regulation .
in the typical case, an over-represented motif did not cover all cluster members but rather a fraction.
the average coverage was 15%, with observations ranging from 1% to 100%.
the limited coverage can be exemplified by the relatively large  cluster of liver-selective genes, which contained  <dig> over-represented motifs at fdr<10%, where the pfm annotation implied hnf <dig>  hnf <dig>  arp <dig>  ppar-γ, and coup-tf as the binding factors .
these different motifs were detected in between  <dig> and  <dig> ortholog pairs, indicating a cluster coverage of less than 20%.
an example of very high coverage was observed in the cluster that contained smooth muscle differentiation markers, in which 6/ <dig>  ortholog pairs were positive for a serum response factor  motif.
we compared our results for hnf- <dig> to a whole-genome experimental screen for targets of this factor in hepatocytes  <cit> .
28%  of the ortholog pairs inside the liver cluster contained experimentally verified hnf- <dig> sites, which indicates a high but not complete coverage.
cross-comparison of experimentally and computationally identified hnf- <dig> targets showed a 69% – 71% agreement .
as a complement to using database motifs, we applied de novo motif elicitation to screen the regulatory sequence of each cluster for over-represented dna motifs that were not present in databases.
a two-step procedure based on the meme algorithm  <cit>  was developed .
the use of a cluster to define a pfm will lead to a bias when testing for over-representation.
a simulation that involved motif elicitation from randomized clusters was used to account for this bias .
<dig> de novo motifs were identified as over-represented in relation to the null simulation .
of these  <dig> motifs,  <dig> matched a described database motif .
9/ <dig> of the matching motifs corresponded to motifs that were over-represented in table  <dig>  the remaining  <dig> motifs may represent novel regulatory elements.
these motifs, and a description of the procedure, are available in the supplementary information.
discussion
the gene battery theory predicts that the core features of differentiated cells are encoded by gene batteries, i.e.
groups of functionally coupled, co-expressed genes that are regulated by similar sets of cis-regulatory elements.
the primary aim of this article was to survey the mouse and human genomes for cases that fall within or close to the idealized gene battery concept.
hierarchical clustering of an extensive compendium of micro-array expression data in mouse and human identified large numbers of co-expressed gene groups.
a grand majority of the larger clusters were significantly enriched for genes that shared one or several go-terms, indicating that co-expressed genes are functionally coupled.
moreover, interactions were significantly more common between proteins that are encoded by clustered pairs compared to randomly chosen pairs.
<dig> clusters, finally, were significantly enriched for genes that shared potential cis-regulatory elements
confirmation and extension of several cases of gene battery-like regulation
the predicted over represented motifs  fell into three principal categories:  predictions in agreement with described gene batteries,  predictions supported by a limited number of observations, where the analysis predicted gene battery-like regulation, and  novel predictions of hypothetical gene batteries.
the findings are summarized with respect to these categories in table  <dig>
for type  predictions the method is useful for proposing new target genes.
as an example, the smooth muscle cluster contains  <dig> validated serum response factor  targets   <cit> .
the method identifies two more genes in this cluster, myl <dig> and lpp as potential srf targets in smooth muscle .
the findings in two clusters that are relevant for lymphocytes  and b-lymphocytes  are further examples of type  predictions.
cluster  <dig> contains genes that are expressed by various lymphocyte populations.
significantly enriched motifs were irf1/irf <dig>  nf-kappab, and isre.
irf <dig> and ir <dig> appear to bind the same sites  <cit> .
irf <dig> knockout mice are impaired in their myelopoiesis  <cit> .
nfkappab is a well-known regulator of inflammation and immune functions  <cit> .
cluster  <dig> contains genes that are specifically expressed in b-lymphocytes.
spib was identified as the most over-represented pfm followed by areb <dig> and oct- <dig>  knockout of spib caused specific defects in b-cell terminally differentiated functions  <cit> .
areb <dig> is a hematopoietic transcriptional repressor  <cit> .
whereas oct- <dig> seems to be dispensable for b-cell development  <cit> , it regulates b lymphocyte genes in combination with specific co-activators  <cit> .
in type  cases, a regulator with some degree of documentation is statistically supported to play a role in gene battery-like regulation.
as an example, rfx family proteins have been shown to regulate a limited number of genes in spermatogenesis, and sox <dig> has been shown to be expressed in the testis  <cit> .
here, we demonstrate that motifs corresponding to these factors are present in 8–11% of testis-selective genes .
29% of kidney-specific genes were covered by hnf- <dig> and hnf- <dig> sites, suggesting that these genes may be general regulators in the kidney.
finally, type  cases represent novel predictions which can be viewed as testable hypotheses of regulatory mechanisms.
a complete list of pfms in the different categories with references to the literature is presented in table  <dig>  examples include a potential functional role for an ikaros-like motif in retina-selective genes , and a pax family factor acting on extracellular matrix genes with strong expression in arteries .
interestingly, a combination of ets family factor motifs, nuclear respiratory factor motifs, and e2f motifs was detected in several clusters of a housekeeping character .
incomplete coverage of cis-regulatory motifs indicates that other mechanisms than co-regulation may contribute to coordinated expression
the gene battery concept predicts that upstream transcription factors coordinate the expression of target genes through binding to similar cis-regulatory elements in the target genes.
in a recent report, alloco and co-workers  <cit>  studied the relationship between co-expression in yeast  and the probability of sharing a transcription factor .
in their experiments, an expression profile correlation of  <dig>  implied a probability of transcription factor sharing of > <dig> .
the incomplete coverage in our analysis may similarly indicate that other mechanisms than co-regulation contribute to the coordinated expression.
importantly, it is technically difficult to accurately quantify the fraction of genes in a cluster that respond to a factor.
the presence of false positive gene cluster members reduces the coverage.
uncertainty in the motif assignment to ortholog pairs may also reduce the coverage.
in one case study, we could use data from an independent experimental screen to confirm that the liver-selective cluster had a limited coverage of hnf- <dig> responders .
further experiments are required to resolve to what extent the lack of coverage reflects alternative regulatory mechanisms or technical limitations.
data limitations
the tissues in the dataset represent samples of limited morphological resolution at a fixed time point.
as a consequence, co-expressed gene groups that are active under specific developmental phases, under specific environmental stresses, or in small and localized anatomical structures are likely to escape detection.
moreover, co-expressed gene groups with a peak in one tissue only are unresolved for single cell types.
low coverage of motifs in such clusters may reflect contaminating genes that derive from another cell type.
the affymetrix technology that were used to generate the expression measurements has been validated and match results obtained with tag sampling for sufficiently abundant genes  <cit> .
low intensity signals did not correlate well, indicating that such genes are less likely to form clusters .
cross hybridization between related genes may, in theory, contribute to correlating expression levels.
we can exclude that clusters primarily form as a result of cross-hybridization, since  members of well-characterized gene families appear in different  clusters  and since  co-expressed genes tended to share go terms or to encode interacting proteins to a higher degree than they tended to be paralogs .
we cannot however exclude subtle effects related to cross-hybridization, since results with cdna microarrays indicate this possibility  <cit> .
transcription factors in the same family often have similar dna binding properties and bind to the same sites on target genes.
this leads to an ambiguity in the interpretation of pfm annotations.
the mef <dig> motif, is for example a receptor for the mammalian mef2a, mef2b, mefc and mefd transcription factors  <cit> .
similarly, the e2f motif is a receptor for a family of  <dig> e2f family proteins that form heterodimers with another family of proteins, the dp proteins, in a way that affects the binding affinity  <cit> .
different e2f proteins act at different stages in the cell cycle  <cit> , and this is not captured by our method.
one additional data limitation is that a fraction of ensembl transcripts may lack sequence in their 5' ends.
the accuracy of our transcription starts is therefore dependent on the quality of the ensembl transcript database.
method considerations and perspectives
our approach is technically related to the toucan and confac tools  <cit> .
important differences and extensions are the introduction of a composite scoring system, a procedure to optimize pfm thresholds, and the use of simulation at the level of the whole clustering to measure significances.
the performance of the method was not affected by the removal of paralogous genes from the clusters , and clustering of such genes is clearly not a significant source of false positive motif prediction.
consequently, we decided to include the paralogs in the final analysis since we believe that co-expression of such genes may be a result of shared cis-regulatory elements.
we further evaluated the effects of using different dna amounts in the analysis.
in principal, it is motivated to include a substantial amount of dna sequence per gene, since mammalian enhancers are frequently located far upstream in relation to the transcription start, or downstream in intronic dna.
however, the benefits of including as much sequence as possible must be balanced against the risk of introducing vast amounts of non-informative sequence into the analysis.
the evaluation clearly favoured using limited amounts of dna, and high stringency phylogenetic footprinting .
the optimal future alternative may be to assemble a small amount of dna from a large area of the genome using algorithms that sort out regulatory dna regions more efficient than phylogenetic footprinting.
the approach presented here needs to be extended to generate more complete models of gene battery regulation.
modelling of cis-regulatory element combinations and the relative position and spacing of elements are examples of such extensions.
the use of pcc as a measure of co-expression may need to be reconsidered if extending the approach to very large datasets, since this measure may be sensitive to noninformative signals in a majority of samples.
further, the expression levels of potential trans-regulators and co-factors  can be introduced.
most likely, the method will benefit from a more accurate identification of regulatory regions
CONCLUSIONS
we screened the mouse and human genomes and transcriptomes for instances of gene battery-like regulation.
comparative clustering was highly predictive of gene function and protein interaction, which indicates that potential gene batteries could be identified this way.
based on a statistical composite score for motifs in ortholog pairs, and a simulation approach to determine significance levels, we found  <dig> instances of statistically supported gene battery-like regulation that were conserved between mouse and human.
these included known cases of gene battery regulation in tissues such as muscle, lymphocytes, erythrocytes, and liver.
a second category of predictions included regulators with some degree of documentation, e.g.
in testis, kidney, and endoplasmatic reticulum.
finally, new candidate gene batteries with statistically enriched cis-regulatory motifs were listed.
the results of this investigation emphasizes the need to study differentiation in terms of larger transcriptional units, and extends the methodology for doing this.

 to understand how gene expression is coordinated to produce hundreds of cell phenotypes from an identical complement of genes is a principal challenge in mammalian genome research. a commonly suggested model for terminal differentiation in metazoans is that the core features of the cellular phenotype are mediated by a set of genes that is regulated as a gene battery, i.e. a set of functionally coupled genes that are activated by similar cis- and trans-acting regulators. although the gene battery is an idealized concept, concrete examples of gene battery-like regulation have been found in for example muscle subtypes, megakaryocytes, the epidermis and lymphocytes. a key step in the elucidation of gene battery-like regulation is to detect and functionally test the cis regulatory elements that mediate the co-regulation. a number of computation-based methods have been proposed to do this. in micro-organisms, computational methods have proven useful to detect modules of co-regulated genes. in mammals, predictive models based on assumed co-regulation at the cis level have been constructed for liver- and skeletal muscle-selective gene regulation, and general tools have been developed for the regulatory analysis of co-expressed genes. the aim of this work is to screen the mouse and human genomes and transcriptomes for instances where sharing of cis-regulatory sequences is statistically coupled to conserved co-expression of genes, i.e. cases that fall within or near the idealized gene battery concept. another aim is to critically investigate technical parameters in order to maximize the sensitivity by which co-regulation of co-expressed genes can be detected. in a screen that covered ~ per cent of all protein-coding genes according to the latest ensembl annotation, we identified co-expressed gene clusters with cases of statistically supported sharing of cis-regulatory sequence motifs. the predictive value of the assignment of transcription factor binding sites was experimentally evaluated on ebf binding sites in a set of b-cell expressed gene clusters. the predicted cases of co-regulation include several previously known prototype examples of tissue specific regulation, but also novel predictions. all data are made available to the research community in the form of an internet resource that may serve as a starting point for further analysis.  the analysis was based on the assumption that homologous genes in mouse and human are equivalent in most aspects of regulation and function. in particular, we assumed that the transcriptional regulation is conserved for orthologous genes. for example, the mouse gene myh and the human gene myh are assumed to share expression pattern and to share important cis-regulatory sequences. below, the term'ortholog pair' will be applied as a two-species equivalent of'gene', for which expression and sequence data where retrieved for both mouse and human. ensembl gene annotations were applied thoughout the analysis. co-expressed gene sets were defined from a compendium of mouse and human expression data and tested for functional coupling previous results by our group and others have shown that statistical analysis of gene expression profiles in a large compendium of expression data can predict targets of differentiation processes, and identify functionally coupled genes. in the following analysis, we specifically focused on a compendium derived from the recently completed novartis expression atlas. these data contain transcription profiles for mouse and human tissues generated by hybridization on customized affymetrix chips, and cover a large fraction of the mouse and human protein-coding genes. sequence annotation of the novartis probes linked the mouse and human data to approximately unique ensembl gene identities in each species. between the two datasets, non-redundant ortholog pairs could be identified by linkage of reciprocal ensembl homology assignments. in later steps, we excluded genes for which regulatory sequence could not be extracted, leaving ortholog pairs for clustering. the final dataset included ~40% of the mouse and human ensembl annotated genes. numbers in the'mouse' and'human' columns signify the number of unique ensembl identifiers in each respective species. numbers in column'both' signify ortholog pairs of ensembl entries. the overlap between the nonredundant sequence database and the nonredundant expression dataset was ortholog pairs. clustering we clustered the mouse/human ortholog pairs based on their expression profiles across the mouse and human tissues. we computed clusters at cut-off levels ranging from pearson's correlation coefficient to. at the lowest applied cut-off, 57% of all ortholog pairs in the data were members of a cluster. the cluster sizes were distributed in a skewed manner, with a predominant formation of small clusters. all analyses hereafter were performed at a pcc = cut-off, which generated clusters with ortholog pairs. this relatively stringent cut-off was chosen to reduce the number of non-relevant genes in the clusters. a higher cut-off did not seem reasonable given the noise-level of the microarray experiments. assessment of functional linkage according to the definition, a gene battery should encode functionally linked proteins. we used gene ontology terms, protein-protein interactions, and manual curation to assess functional linkage within our expression clusters. first, we investigated the relationship between expression profile similarity of two ortholog pairs and their relative probability to share a functional annotation term or to encode interacting proteins. there was a consistent correlation between co-expression and the relative probability for two genes to share a go-term or to encode interacting proteins. second, we studied the statistical over-representation of go terms and interacting proteins inside the clusters. 30/ clusters with ten or more ortholog pairs contained at least one over-represented go term. the proportion of small clusters with over-represented go terms was lower, which reflects a lack of statistical power in small clusters. genes encoding interacting protein pairs were also over-represented inside clusters. cases of protein-protein interaction between two genes in the same cluster were found. in contrast, simulations on permuted data revealed a median of interactions. the bind data contained only interactions that could be mapped to the dataset, which explains the seemingly low number of interactions. finally, the clusters were annotated by manual curation. six examples are shown in figures and clusters with predicted regulators are listed in table for a full overview, see the web supplement. several clusters clearly represent gene sets that mediate specialized features of different cell types, including smooth muscle specific genes, b lymphocyte-specific genes, and genes selectively expressed in the testis. further, we detected clusters that were related to cellular processes or organelles, including endoplasmatic reticulum, protein synthesis, and modulators of transcriptional regulation. a grand majority of the clusters were defined by peaks at different amplitudes in several tissues. as an example, the endoplasmatic reticulum cluster was defined by a highly variant profile with strong expression in, for example, exocrine glands. generally, the cluster profiles were conserved between species, in the sense that clusters were defined by peaks in the same organs. this effect was more pronounced for clusters with expression in a single organ, such as the testis. over-represented motifs arranged by cluster number. fdr column: false discovery rate. motifs are are shown both by their numerical identifiers and by their annotation. in cases where a pfm is a composite based on more than one source, the components are given separated by commas. the data where generated from the pcc = clustering, kb sequence database, at 90% phylogenetic conservation. in combination, the go term enrichment, the protein-protein interaction, and the manual curation convincingly show that clustered genes are functionally linked. regulatory dna and descriptions of transcription factor binding sites were extracted and pre-processed in the next part of the analysis, individual ortholog pairs were scored for transcription factor binding sites. binding motifs were represented in the form of position frequency matrices. based on a fixed amount of upstream dna sequence in each ortholog pair, a statistical score was computed to predict the potential for a site in the sequence to bind the factor corresponding to a pfm. we extracted upstream + intronic sequence from the ensembl database, in amounts of or kb per gene. a filter was applied that removed ortholog pairs for which the transcription start differed between the two species. filtering was successful for unique ortholog pairs. in a subsequent filter, dna sequence that was not conserved between mouse and human was removed, so-called phylogenetic footprinting. phylogenetic footprinting was applied at different stringency, to allow the following optimization of the protocol. furthermore, all exon sequence was removed from the analysis. finally, the sequences were matched with the expression data based on annotation, the overlap being ortholog pairs. vertebrate pfms were downloaded from the transfac and jaspar databases. since the databases appeared to contain redundant or equivalent entries, highly similar pfms were grouped and merged using single linkage hierarchical clustering and a pfm distance measure defined in, which reduced the number of pfms from to. this step reduced the redundancy, but did not merge all identically annotated pfms pfms in table cluster 40). design of a predictive scoring system for transcription factor binding after the retrieval and preprocessing of both sequences and pfms, all individual sequences were tested for pfm matches using the mast software, a software for identifying single or multiple motifs in sequences. mast was set to compute one p-value for each pfm with respect to each sequence. based on the p-values obtained from the mast software, a composite score was defined as the product of the p-value in the mouse and human sequences of an ortholog pair. a composite score close to indicates that both the mouse and human promoter sequence in the ortholog pair contains sequence elements that are in very good agreement with a certain pfm. to address the biological validity of the mast composite scores within the context of a set of co-expressed genes, we screened ortholog pairs present in b-cell expressed clusters for individual ebf sites. in all, individual ebf sites in different ortholog pairs were detected. to test the functionality of these sites, we examined the potential of basepair duplex oligo-nucleotides spanning the presumed sites to compete for protein binding in electrophoretic mobility shift assays. ebf binding capacity was assayed using nuclear extracts from the pre-b cell line 40-ei, a labelled mouse mb- promoter high affinity ebf site, and competitor oligo-nucleotides covering the new potential binding sites. in the absence of competitor oligo-nucleotide, a prominent dna/protein complex could be detected whereas this complex was undetectable after the addition of a fold molar excess of the unlabeled binding site. the identity of the protein and the specificity of the binding were verified by competition with a point mutated mb-1/ebf site and by the inclusion of an ebf specific antibody into the reaction mixture. the point mutated ebf binding site was unable to abolish complex formation even in a 1000-fold molar excess, indicating that we detect specific protein dna interactions with this experimental set up. out of the new binding sites competed for complex formation when added in a 300- or 1000-fold molar excess, and thus have the ability to bind ebf in vitro. we conclude that the large majority of binding sites were able to functionally interact with the predicted protein, and that the composite score in principle detected factor binding. it should however be emphasized that the quality of predictions is dependent on the quality of the binding site descriptions, and the result does not necessarily imply that other predicted factors bind. a statistical procedure was used to detect enriched motifs in the clusters to test whether the identified clusters represent potential gene batteries, i.e. contain shared cis-regulatory elements, we designed a procedure to detect significant over-representation of orthologs that match a pfm inside a cluster. the procedure is based on a modification of fisher's exact test, which tests for dependency between two events. we introduced a procedure to optimize the composite score thresholds for individual pfms. in brief, we selected the threshold that gave the lowest fisher test p-value in any one cluster. this was based on the assumption that non-random distribution of detections over clusters reflects biological function, as has been proposed in. the fisher test p-values with the optimized thresholds are hereafter termed p-scores. the tests of multiple detection thresholds, multiple clusters and pfms led to a need to compensate for mass testing. this was done by estimating false discovery rates based on simulations on randomized data. in brief, we compared the outcome when using permuted and when using observed data at different p-scores, and defined the fdr as the ratio between the two. this procedure allowed us to choose a significance threshold with a controlled expected number of spurious detections. simulations were repeated times. since optimization of detection thresholds were repeated in each simulation round, no bias in disfavour of the control case was introduced. the detection of over-represented pfms was affected by dna amount and masking, but not affected by gene paralogy using the described algorithm, we analyzed the number of times a motif was over-represented when using different amounts of dna sequence per gene and different stringency in the phylogenetic footprinting, >60%, >70%, >80% or >90% identity). the analysis was performed at 10% false discovery rate and with the clustering obtained at pcc = clustering cut-off. at all dna amounts, higher stringency phylogenetic footprinting appeared to be beneficial, and between dna amounts, kb and kb compared favourably over kb. we conclude that optimal results are obtained when using a limited amount of sequence per gene. a potential confounding factor in the analysis is that similarity in upstream sequence may be attributable to factors other than shared cis-regulatory elements. the most important such factor is likely to be gene paralogy, since over-represented motifs might simply represent matches to non-functional sequence in a co-expressed gene family. in total, cases of pairwise paralogy were detected in the dataset of ortholog pairs using blast. all genes with transcripts matches at a blast e-value less than were defined as potential paralogs. of the ortholog pairs that clustered, genes had a paralog inside the same cluster. this was clearly an over-representation since filling clusters with random genes from the dataset produced an average of genes with a paralog inside a cluster. in a side-by-side comparison, we therefore analyzed the effect of removing paralogs from a cluster before testing for over-represented motifs, as opposed to not removing paralogs but the same number of randomly selected genes. the results showed no increase in the number of over-represented motifs when allowing paralogs in the same cluster. over-represented motifs were observed in clusters the final analysis was performed with kb sequence and 80% phylogenetic footprinting. again, the clustering generated at pcc = was chosen for analysis. the over-representation algorithm was run in the above cases, using simulations to estimate false discovery rate thresholds. over-represented motifs with fdr:s less than % and 10% were recorded. key features of these results are presented in tables and and the complete results are available by web browser. criteria for inclusion in the different groups: i. existing evidence of co-regulation of cluster member by the factor, and/or a mutant mouse phenotype that affects differentiation of the cell type. ii. evidence of the factor acting on a limited number of genes in the tissue. evidence of the transcription factor being expressed in the tissue. iii. no or limited evidence in the literature. * = detected at 20% fdr. a predictive screen performed with these settings associated motifs to co-expression in a total of clusters at fdr thresholds % and 10%. there was a clear tendency for motifs to be detected in most large clusters, and for smaller clusters to lack over-represented motifs. this can be directly explained by the fact that even a high degree of motif presence in a small cluster can be attributed to spurious detections, and that the fisher test accounts for that. the over-represented motifs were validated against the literature, and fell into three main categories: 1) previously described cases of gene-battery-like regulation, 2)previously unreported cases of gene battery-like regulation with some support in a limited number of genes, and 3) hypothetical cases of gene battery regulation. in the typical case, an over-represented motif did not cover all cluster members but rather a fraction. the average coverage was 15%, with observations ranging from 1% to 100%. the limited coverage can be exemplified by the relatively large cluster of liver-selective genes, which contained over-represented motifs at fdr<10%, where the pfm annotation implied hnf hnf arp ppar-γ, and coup-tf as the binding factors. these different motifs were detected in between and ortholog pairs, indicating a cluster coverage of less than 20%. an example of very high coverage was observed in the cluster that contained smooth muscle differentiation markers, in which 6/ ortholog pairs were positive for a serum response factor motif. we compared our results for hnf- to a whole-genome experimental screen for targets of this factor in hepatocytes. 28% of the ortholog pairs inside the liver cluster contained experimentally verified hnf- sites, which indicates a high but not complete coverage. cross-comparison of experimentally and computationally identified hnf- targets showed a 69% – 71% agreement. as a complement to using database motifs, we applied de novo motif elicitation to screen the regulatory sequence of each cluster for over-represented dna motifs that were not present in databases. a two-step procedure based on the meme algorithm was developed. the use of a cluster to define a pfm will lead to a bias when testing for over-representation. a simulation that involved motif elicitation from randomized clusters was used to account for this bias. de novo motifs were identified as over-represented in relation to the null simulation. of these motifs, matched a described database motif. 9/ of the matching motifs corresponded to motifs that were over-represented in table the remaining motifs may represent novel regulatory elements. these motifs, and a description of the procedure, are available in the supplementary information. discussion the gene battery theory predicts that the core features of differentiated cells are encoded by gene batteries, i.e. groups of functionally coupled, co-expressed genes that are regulated by similar sets of cis-regulatory elements. the primary aim of this article was to survey the mouse and human genomes for cases that fall within or close to the idealized gene battery concept. hierarchical clustering of an extensive compendium of micro-array expression data in mouse and human identified large numbers of co-expressed gene groups. a grand majority of the larger clusters were significantly enriched for genes that shared one or several go-terms, indicating that co-expressed genes are functionally coupled. moreover, interactions were significantly more common between proteins that are encoded by clustered pairs compared to randomly chosen pairs. clusters, finally, were significantly enriched for genes that shared potential cis-regulatory elements confirmation and extension of several cases of gene battery-like regulation the predicted over represented motifs fell into three principal categories: predictions in agreement with described gene batteries, predictions supported by a limited number of observations, where the analysis predicted gene battery-like regulation, and novel predictions of hypothetical gene batteries. the findings are summarized with respect to these categories in table for type predictions the method is useful for proposing new target genes. as an example, the smooth muscle cluster contains validated serum response factor targets. the method identifies two more genes in this cluster, myl and lpp as potential srf targets in smooth muscle. the findings in two clusters that are relevant for lymphocytes and b-lymphocytes are further examples of type predictions. cluster contains genes that are expressed by various lymphocyte populations. significantly enriched motifs were irf1/irf nf-kappab, and isre. irf and ir appear to bind the same sites. irf knockout mice are impaired in their myelopoiesis. nfkappab is a well-known regulator of inflammation and immune functions. cluster contains genes that are specifically expressed in b-lymphocytes. spib was identified as the most over-represented pfm followed by areb and oct- knockout of spib caused specific defects in b-cell terminally differentiated functions. areb is a hematopoietic transcriptional repressor. whereas oct- seems to be dispensable for b-cell development, it regulates b lymphocyte genes in combination with specific co-activators. in type cases, a regulator with some degree of documentation is statistically supported to play a role in gene battery-like regulation. as an example, rfx family proteins have been shown to regulate a limited number of genes in spermatogenesis, and sox has been shown to be expressed in the testis. here, we demonstrate that motifs corresponding to these factors are present in 8–11% of testis-selective genes. 29% of kidney-specific genes were covered by hnf- and hnf- sites, suggesting that these genes may be general regulators in the kidney. finally, type cases represent novel predictions which can be viewed as testable hypotheses of regulatory mechanisms. a complete list of pfms in the different categories with references to the literature is presented in table examples include a potential functional role for an ikaros-like motif in retina-selective genes, and a pax family factor acting on extracellular matrix genes with strong expression in arteries. interestingly, a combination of ets family factor motifs, nuclear respiratory factor motifs, and e2f motifs was detected in several clusters of a housekeeping character. incomplete coverage of cis-regulatory motifs indicates that other mechanisms than co-regulation may contribute to coordinated expression the gene battery concept predicts that upstream transcription factors coordinate the expression of target genes through binding to similar cis-regulatory elements in the target genes. in a recent report, alloco and co-workers studied the relationship between co-expression in yeast and the probability of sharing a transcription factor. in their experiments, an expression profile correlation of implied a probability of transcription factor sharing of >. the incomplete coverage in our analysis may similarly indicate that other mechanisms than co-regulation contribute to the coordinated expression. importantly, it is technically difficult to accurately quantify the fraction of genes in a cluster that respond to a factor. the presence of false positive gene cluster members reduces the coverage. uncertainty in the motif assignment to ortholog pairs may also reduce the coverage. in one case study, we could use data from an independent experimental screen to confirm that the liver-selective cluster had a limited coverage of hnf- responders. further experiments are required to resolve to what extent the lack of coverage reflects alternative regulatory mechanisms or technical limitations. data limitations the tissues in the dataset represent samples of limited morphological resolution at a fixed time point. as a consequence, co-expressed gene groups that are active under specific developmental phases, under specific environmental stresses, or in small and localized anatomical structures are likely to escape detection. moreover, co-expressed gene groups with a peak in one tissue only are unresolved for single cell types. low coverage of motifs in such clusters may reflect contaminating genes that derive from another cell type. the affymetrix technology that were used to generate the expression measurements has been validated and match results obtained with tag sampling for sufficiently abundant genes. low intensity signals did not correlate well, indicating that such genes are less likely to form clusters. cross hybridization between related genes may, in theory, contribute to correlating expression levels. we can exclude that clusters primarily form as a result of cross-hybridization, since members of well-characterized gene families appear in different clusters and since co-expressed genes tended to share go terms or to encode interacting proteins to a higher degree than they tended to be paralogs. we cannot however exclude subtle effects related to cross-hybridization, since results with cdna microarrays indicate this possibility. transcription factors in the same family often have similar dna binding properties and bind to the same sites on target genes. this leads to an ambiguity in the interpretation of pfm annotations. the mef motif, is for example a receptor for the mammalian mef2a, mef2b, mefc and mefd transcription factors. similarly, the e2f motif is a receptor for a family of e2f family proteins that form heterodimers with another family of proteins, the dp proteins, in a way that affects the binding affinity. different e2f proteins act at different stages in the cell cycle, and this is not captured by our method. one additional data limitation is that a fraction of ensembl transcripts may lack sequence in their 5' ends. the accuracy of our transcription starts is therefore dependent on the quality of the ensembl transcript database. method considerations and perspectives our approach is technically related to the toucan and confac tools. important differences and extensions are the introduction of a composite scoring system, a procedure to optimize pfm thresholds, and the use of simulation at the level of the whole clustering to measure significances. the performance of the method was not affected by the removal of paralogous genes from the clusters, and clustering of such genes is clearly not a significant source of false positive motif prediction. consequently, we decided to include the paralogs in the final analysis since we believe that co-expression of such genes may be a result of shared cis-regulatory elements. we further evaluated the effects of using different dna amounts in the analysis. in principal, it is motivated to include a substantial amount of dna sequence per gene, since mammalian enhancers are frequently located far upstream in relation to the transcription start, or downstream in intronic dna. however, the benefits of including as much sequence as possible must be balanced against the risk of introducing vast amounts of non-informative sequence into the analysis. the evaluation clearly favoured using limited amounts of dna, and high stringency phylogenetic footprinting. the optimal future alternative may be to assemble a small amount of dna from a large area of the genome using algorithms that sort out regulatory dna regions more efficient than phylogenetic footprinting. the approach presented here needs to be extended to generate more complete models of gene battery regulation. modelling of cis-regulatory element combinations and the relative position and spacing of elements are examples of such extensions. the use of pcc as a measure of co-expression may need to be reconsidered if extending the approach to very large datasets, since this measure may be sensitive to noninformative signals in a majority of samples. further, the expression levels of potential trans-regulators and co-factors can be introduced. most likely, the method will benefit from a more accurate identification of regulatory regions  we screened the mouse and human genomes and transcriptomes for instances of gene battery-like regulation. comparative clustering was highly predictive of gene function and protein interaction, which indicates that potential gene batteries could be identified this way. based on a statistical composite score for motifs in ortholog pairs, and a simulation approach to determine significance levels, we found instances of statistically supported gene battery-like regulation that were conserved between mouse and human. these included known cases of gene battery regulation in tissues such as muscle, lymphocytes, erythrocytes, and liver. a second category of predictions included regulators with some degree of documentation, e.g. in testis, kidney, and endoplasmatic reticulum. finally, new candidate gene batteries with statistically enriched cis-regulatory motifs were listed. the results of this investigation emphasizes the need to study differentiation in terms of larger transcriptional units, and extends the methodology for doing this. 
 predicted cases of over-represented transcription factor binding motifs were validated against the literature and fell into three categories: previously described cases of gene battery-like regulation, previously unreported cases of gene battery-like regulation with some support in a limited number of genes, and predicted cases that currently lack experimental support. the novel predictions include for example sox and rfx transcription factor binding sites that were detected in ~10% of all testis specific genes, and hnf- and binding sites that were detected in ~30% of all kidney specific genes respectively. we developed a predictive procedure to screen the mouse and human genomes and transcriptomes for cases of gene-battery-like regulation. the expression of gene batteries, genomic units of functionally linked genes which are activated by similar sets of cis- and trans-acting regulators, has been proposed as a major determinant of cell specialization in metazoans. in a screen that covered ~ per cent of all annotated protein-coding genes, we identified co-expressed gene clusters with statistically supported sharing of cis-regulatory sequence elements. co-expressed gene clusters were enriched for a total of shared cis-regulatory sequence elements. the results are publicly available at. critical technical parameters were evaluated, and the results and the methods provide a valuable resource for future experimental design.
 the expression of gene batteries, genomic units of functionally linked genes which are activated by similar sets of cis- and trans-acting regulators, has been proposed as a major determinant of cell specialization in metazoans. we developed a predictive procedure to screen the mouse and human genomes and transcriptomes for cases of gene-battery-like regulation.  in a screen that covered ~ per cent of all annotated protein-coding genes, we identified co-expressed gene clusters with statistically supported sharing of cis-regulatory sequence elements. predicted cases of over-represented transcription factor binding motifs were validated against the literature and fell into three categories: previously described cases of gene battery-like regulation, previously unreported cases of gene battery-like regulation with some support in a limited number of genes, and predicted cases that currently lack experimental support. the novel predictions include for example sox and rfx transcription factor binding sites that were detected in ~10% of all testis specific genes, and hnf- and binding sites that were detected in ~30% of all kidney specific genes respectively. the results are publicly available at.  co-expressed gene clusters were enriched for a total of shared cis-regulatory sequence elements. a majority of these predictions represent novel cases of potential co-regulation of functionally coupled proteins. critical technical parameters were evaluated, and the results and the methods provide a valuable resource for future experimental design. 