BACKGROUND
variation in gene expression provides a quantifiable trait that has been employed to classify breast tumors  <cit> .
however it has long been known that the gene sets identified from independent laboratories fail to provide a unified set of genes thereby casting doubt on the biological implications of these profiles  <cit> .
despite these differences, two prognostic tests have recently been approved in the united states for clinical management of disease  <cit> .
from a diagnostic perspective, developing a unified gene profile that predicts both risk of recurrence and therapeutic response in diverse disease subtypes would be clinically useful.
these gene sets could also provide an understanding of the mechanistic basis of malignancy.
meta-analysis has been used as a formal summarization method in the clinical cancer literature for many years  <cit> .
recently, some groups have applied meta-analysis to gene expression microarrays  <cit> .
meta-analysis refers to a broad class of models used for summarizing and synthesizing studies to estimate their overall effect.
rhoades, et al was among the first to demonstrate the usefulness of meta-analytic procedures on microarray data in prostate cancer  <cit> .
since then, there have been many contributions to the oncology literature by applying meta-analysis to microarrays, including breast cancer  <cit> .
one of the central goals in gene expression experiments is to identify the common regulatory themes and cis-elements responsible for the observed patterns of gene expression.
this has been most successfully performed for the yeast saccharomyces cerevisiae where new regulatory genes have been suggested  <cit> .
however, metazoan expression patterns tend to be more complicated.
one approach has been to combine expression data of orthologous genes from diverse organisms to build co-expression networks  <cit> .
in drosophilia gene networks have been proposed based upon the co-localization of tfbs with cis-regulatory modules   <cit> .
the availability of both mammalian and lower metazoan complete genome assemblies affords one the opportunity to identify phylogenetically conserved motifs in the array candidates.
in addition to known tfbs, these phylogenetic motifs may identify important new cis-acting signals that modulate transcription  or transcript stability  and may be key elements in the observed expression patterns.
a systematic comparison of both known and phylogenetic cis-elements between two sets of differentially expressed genes can serve to implicate these elements as common modulators in the observed gene expression patterns.
our method incorporates a meta-analysis model to rank genes into groups of over- and under-expressed gene sets, based upon their relative importance between independent array studies.
our analyses of gene expression patterns in er+ and er- breast tumors were performed across different array platforms on a diverse spectrum of patients.
the two sets of genes showing the most disparate expression patterns between er+ and er- tumors provided an entry point with which to explore the possibility that specific sequence elements may be disproportionately represented in these two groups.
we utilized known motifs in conjunction with comparative genomic resources to search for enriched dna elements in both the proximal promoter and 3'utr regions of these genes.
our findings suggest that the differential gene expression in er+ vs. er- tumors may, in some cases, be mediated by specific sequence elements in either the promoter or 3'utr intervals.
the motif distribution profiles between our gene sets identified both known and phylogenetically conserved elements that may play a role in these genes' co-expression.
RESULTS
forty-six percent of unique probes among the studies mapped many-to-one to unique unigene ids.
the mean and median numbers of probes per unigene ids were  <dig>  and  <dig>  respectively.
when we merged the  <dig> studies in table  <dig> for the meta-analysis data set, we retained the expression values for all probe combinations in all studies and this resulted in a multiplicative set of records in the database.
approximately 12% of the unique ests in the oncomine database  did not correspond to a unique unigene id.
these were dropped from the analysis data sets.
dfs, disease free survival; pr, progesterone status; ln, lymph node status; os, overall survival
we focused our subsequent analyses on a select set of genes by taking medians across each unigene id's s/sd statistics.
a scatter plot of the s/n  versus abs/sd statistics  appears in figure  <dig>  the distribution of the s/n values was bell-shaped with heavy tails.
our criteria for selecting genes were to take the most extreme 1% and 5% values in both tails.
we found it instructive to consider the ratio s/sd on the y-axis of figure  <dig>  where sd is the standard deviation of the  addends of s. large values of this ratio indicate those genes with consistently significant p-values across all of the studies that we considered.
the number of unigene ids with s/n scores in the top 1% and 5%  were  <dig> and  <dig>  respectively.
the mean numbers of studies for genes present in the top 1% and top 5% classes were  <dig>  and  <dig>  respectively.
our choice of reporting both top 1% and the top 5% for further screening was for crude management of false positives from bias correlated with each gene's relative ranking.
many of the genes present in our top 1% upregulated list identified in our meta-analyses have previously been identified as overexpressed in er+ breast tumors, most notably the two transcription factors esr <dig> and gata <dig>  our gene lists appear in additional file  <dig>
we next compared our top 1% and 5% upregulated gene lists in er+ and er- tumors to those prognostic genes utilized in the 70-gene signature associated with mammaprint®  <cit>  along with the 16-gene signature with the rt-pcr based oncotype dx®  <cit>  tests.
although the array data defining the 70-gene profile was one of  <dig> input datasets for our meta-analysis and the 16-gene signature datasets utilized two expression datasets were also employed for the analysis, we did not observe complete overlap in the genes identified.
for the 70-gene signature our top 1% dataset identified an overlap of one and four genes respectively that were upregulated in er+ tumors versus er- tumors.
only  <dig> and  <dig> genes overlapped in the top 5% dataset, respectively.
alternatively, for the 16-gene signature, one and two genes, respectfully, from the top 1% gene sets were overexpressed in er+ versus er- tumors from our meta-analysis, while  <dig> and  <dig> genes, respectfully, overlapped in the top 5% list.
differences in probes, arrays, and studies used in the meta-analyses may explain some of the differences between our gene lists and the gene lists from the two diagnostic tools.
additionally, we compared our gene lists to a previously identified universal profile that uses  <dig> genes overexpressed in a diverse spectrum undifferentiated cancers to predict neoplastic transformation  <cit> .
strikingly we observed only genes overexpressed from er- tumors to overlap with this  <dig> gene signature.
four genes  overlapped in our top 1% list while  <dig> genes overlapped from our top 5% list.
these findings are in accord with the observation that er- tumors are more highly undifferentiated than er+ tumors and demonstrate more metastatic potential clinically  <cit> .
collectively these  <dig> genes may identify a set of candidate genes predictive of metastatic potential in er- breast tumors.
ingenuity pathway analyses
we considered the relationship of our top 1% genes in the er+ and er- groups using ingenuity pathway analyses  <cit> .
our objective in using ingenuity was to characterize the functional role of our selected genes.
ipa isolated genes for which it had documented associations, and created a series of networks based on the published literature.
we were able to map  <dig> of the  <dig> genes comprising the sum of the 1% upregulated and 1% downregulated gene sets.
from these networks, ipa queried its database of biological functions and scored each gene cluster with a p-value calculation.
table  <dig> shows the most common functions found among our most differentially-expressed genes.
notably our top 1% genes upregulated in er- tumors contained  <dig> genes showing association to cancer whereas only  <dig> of the genes upregulated in er+ tumors were cancer-associated.
promoter motif comparisons in dysregulated genes
we tested the hypothesis that there was a significant difference in the occurrence of each motif between our two classes of genes  using a fisher's exact test.
we adjusted for multiple testing by applying the benjamini-hochberg p-value correction  <cit> .
we counted the number of genes in each class which were overexpressed in er+ tumors and contained a copy of each phylogenetic motif, and compared those to the number of genes overexpressed in er- tumors.
for genes harboring multiple copies of a motif we counted these elements as a single motif event.
we independently performed tests for both the top 1% and 5% of our genes.
our initial query sets consisted of  <dig> condensed transfac motifs and a second analysis comprised  <dig> phylogenetically conserved mammalian promoter motifs as previously defined  <cit> .
sixty-nine of the phylogenetic motifs map to known tfbs defined in the transfac db v <dig>  while  <dig> represented novel phylogenetically conserved elements.
we first examined whether any of  <dig> known tfbs were disproportionately represented in our er+ and er- gene sets.
abbreviated results appear in table  <dig>  while numerous motifs showed significance by fisher's exact testing  only  <dig> survived multiple testing correction.
the first motif ktwgttt, a binding site for the sry <dig> transcription factor, was over-represented in the top 5% of er+ upregulated genes in the noncoding strand.
for er+ overexpressed genes  <dig> of  <dig> genes contained the site while  <dig> of  <dig> er- overexpressed genes contained the site .
the second site, abwcaggtrnr, a binding site for areb <dig> , was over-represented in the top 1% of er+ upregulated genes when both coding and noncoding strands were surveyed  and contains an embedded e-box motif.
twenty-five genes bore tcf <dig> sites in either strand amongst  <dig> er+ upregulated genes while only  <dig> genes contained the site amongst  <dig> er- upregulated genes.
the presence of tcf <dig> sites in nearly four times as many er+ upregulated genes versus er- upregulated genes may be an indirect mechanism for gene activation in er+ breast tumors.
tcf <dig> has been shown to be induced by estrogen which in turn activates a cascade of downstream genes  <cit> .
additionally, the transcriptional repression of e-cadherin by tcf <dig> has been shown to lead to loss of the epithelial phenotype suggesting a role for this tf in late-stage carcinogenesis  <cit> .
we note that although e-cadherin was not identified in our meta-analysis,  <dig> related genes, cdh <dig> and pcdh <dig>  both of which lie in the top 5% of er- overexpressed genes, may be responsive to repression by tcf <dig>  the over-representation of tcf <dig> binding sites in both strands of our top 1% genes er+ overexpressed tumors suggests that tcf <dig> may act as a transcriptional activator for these genes yet act as a transcriptional repressor in er- overexpressed genes.
in addition to known sites, we sought to identify potential new regulatory motifs by examining the coding and noncoding strands with  <dig> previously identified phylogenetic motifs in the top 1% and 5% of our s+ vs. s- genes  <cit> .
eleven of these motifs represented palindromic sequences and were scanned in only the coding strand when both strands were analyzed.
again, while numerous motifs showed significance by fisher's exact testing  only  <dig> survived multiple testing correction.
abbreviated results appear in table  <dig>  a single motif  showed a significant difference between the er+ upregulated genes versus the er- upregulated genes when the non-coding strand was examined in our top 1% gene list.
nineteen of  <dig> er+ overexpressed genes contained at least  <dig> copy of the motif while only  <dig> of  <dig> genes contained the motif in the er- overexpressed genes .
this phylogenetic motif does not map to any known tfbs and represents a new target for exploration.
*center of motif relative to transcriptional start site
analysis of the 3'utr
we next screened for regulatory elements in the 3'utr of our genes sets.
less is known about functional motifs in 3'utrs than about functional motifs in promoter regions, but evolutionary conserved motifs in 3'utrs may, as in promoter regions, indicate regulatory sites.
we therefore used a previously identified set of evolutionary conserved 3'utr motifs  <cit> .
although the function of half of these motifs is unknown, the remaining half has a/t rich elements believed to be involved in mrna stability or represent likely micrornas binding sites.
we used the same refseq id's to harvest the annotated 3'utrs of our gene sets as described in the methods.
surprisingly, we observed a significant difference in the median 3'utr lengths between our gene sets .
the top 1% genes overexpressed in er+ tumors contained a median 3'utr length of  <dig>  kb, while genes overexpressed in er- tumors contained a median 3'utr length of  <dig>  kb.
a similar trend was observed when we examined the top 5% of genes sets.
in this set, er+ upregulated genes had a median utr length of  <dig>  kb while the er- genes had a length of  <dig>  kb.
microrna target genes have longer 3'utrs, whereas anti-targets have shorter 3'utrs  <cit> .
thus, the difference in 3'utr length suggests a difference in mirna targeting prevalence between the er+ and er- genes.
the most significant evolutionary conserved motif in the top 1% and top 5% genes  correspond to a potential mirna target site; yactgccr and wgcctta have seed complementarity to mir-34/mir- <dig> and mir- <dig>  the mirna seed region – nucleotides 2– <dig> from the 5' end – is the most important factor for mirna target site recognition  <cit> .
fisher's exact tests on the mirna seed site occurrence counts, corrected for multiple testing, seemingly confirm that the er+ genes are preferentially regulated by mirnas, as all the significant seeds are overrepresented in the er+ upregulated genes .
there is, however, a potential problem with using the fisher's exact test for the 3'utr sets.
if motif occurrences were random, we would expect the er+ genes to have more motif occurrences than the er- genes have, as the er+ genes have longer 3'utrs.
thus, to determine whether there is a significant difference in mirna regulation between the er+ and er- genes, we had to address whether the occurrences of mirna seed sites in the two sets were significantly different from what we would expect by chance.
we therefore ran a set of randomization experiments where we compared the observed number of seed site occurrences in the er+ and er- genes' 3'utrs with those in random gene sets that had 3'utr lengths similar to the er+ and er- 3'utrs .
we found that all of the seeds identified by significant fisher's exact tests do occur significantly more frequently in the er+ 3'utrs than in 3'utrs from random gene sets.
moreover, these seeds also occur significantly less frequently in the er- 3'utrs than in random gene sets.
thus, it seems that whereas several mirnas may coordinately regulate some of the er+ genes, some of the er- genes may collectively avoid being regulated by the same mirnas.
previous studies have identified several mirnas that are aberrantly expressed in breast cancers  <cit> .
together the aberrantly expressed mirnas in these studies mapped to  <dig> unique  <dig> mer seed sequences of which three were among our ten most significant  <dig> mer motifs.
the three corresponding mirnas  are all overexpressed in breast cancers.
none of the ten most significant  <dig> mer motifs are from mirnas reported to be differentially expressed in er+ and er- tumors  <cit> ; the most significant  <dig> mer is ranked 25th , is significantly more abundant in er+ genes than expected by random, and is from mir- <dig>  which is downregulated in er+ tumors.
discussion
micrornas are small  noncoding rnas that recognize complementary target sequences in mrnas and prompt either translational repression or rna degradation.
micrornas play important roles in cancer.
iorio et al., for example, recently revealed that deregulation of multiple mirnas can be correlated to pathogenic features such as estrogen or progesterone receptor status and tumor stage for breast cancers  <cit> .
in addition, shorter postoperative survival times for patients with lung tumors can be predicted by measuring mirna let- <dig>  <cit> .
thus mirnas can be used both as classifiers of breast tumor type and as predictors of survival of lung cancer patients.
micrornas preferentially target 3'utrs that have short sequences with perfect complementarity to nucleotides 2– <dig>  or 2– <dig>  in the mirna's 5' region – the seed region  <cit> .
as mirna regulation may explain gene co-expression, we therefore included the  <dig> mer and  <dig> mer seed sequences for all human mirna sequences known at the time of the study.
we note that not all known human mirnas are highly evolutionary conserved and these seed sequences therefore supplement the mirna-related evolutionary conserved motifs.
since we identified sets of genes that demonstrated differential expression between er+ and er- tumors, we reasoned that some of these genes may contain common cis-regulatory motifs contributing to their co-regulation.
we would predict that these sites may, in some cases, be disproportionately represented between genes upregulated in er+ tumors versus genes upregulated in er- tumors perhaps allowing one to identify genes sharing common regulatory pathway.
computational tools exist to identify tfbs based upon over-representation of conserved motifs in datasets  <cit> .
other approaches aim to identify transcription factors  which bind to tfbs based on the relatedness of expression profiles between the tf and the target genes they are postulated to regulate  <cit> .
a combined approach utilizing expression measurements of tissue-specific gene sets in conjunction with orthologous tfs from humans and mouse provides for enhanced accuracy in predicting bone fide cis-regulatory elements  <cit> .
for the most part these searches are guided by biologically confirmed tfbs interactions identified in the transfac database  <cit> ; however, this approach may fail to identify motifs that may be evolutionarily conserved amongst mammals.
in addition to known sites that remained significant after multiple testing correction, many additional sites, and their associated transcription factors, warrant comment.
a second important tfbs, ctttga, the binding site for lymphoid enhancer-binding factor  <dig> , in the top 1% coding strand er+ overexpressed genes, failed rigorous multiple testing where  <dig> of  <dig> genes contained ≥  <dig> site versus  <dig> of  <dig> genes in er- gene set in table  <dig>  nonetheless there is strong biological evidence supporting the role of lef <dig> in tumorogenesis.
the lef <dig> binding site ctttga is one of the primary binding sites in the wnt signaling pathway which regulates cell-cell adhesion and many morphogenetic events during mammary development and possibly cancer  <cit>  binding of wnt proteins with frizzled protein prevents degradation of β-catenin, which subsequently translocates to the nucleus and binds transcription factors of the tcf/lef family .
several tumors are known to have an altered β-catenin signaling pathway including colorectal and lymphoblastic tumors  <cit> .
mutations in the wnt pathway genes can result in β-catenin stabilization and activation of lef/tcf-induced transcription.
recent studies have demonstrated sebaceous tumors harboring lef <dig> mutations interfere with β-catenin-binding domain of lef <dig> and transcriptional activation  <cit> .
common human carcinomas also carry mutations in the β-catenin-binding domain of lef  <cit> .
our data suggest that mutations  in lef <dig> or tcf <dig> binding sites in genes that inactivate wnt signaling could contribute to breast tumorogenesis.
we did not find the estrogen receptor binding site  over-enriched in any our analyses.
this is not surprising as our survey was confined to the immediate  <dig> kb promoter region.
we point out that estrogen may be playing an indirect role on genes in er+ overexpressing tumors via the activation of tf such as tcf <dig> which in turn activate downstream targets.
additionally, it is possible that differences in er binding sites do exist between our gene sets but these sites may reside at distances much further upstream.
recent reports indicate that only two-thirds of er tfbs can be localized to the proximal promoter region of rna polymerase ii genes  <cit> .
we also note that the e2f binding site  consistently ranked amongst the top  <dig> motifs  identified when screening the non-coding strand.
in the non-coding strand of the top 1% gene sets, more e2f sites were observed in genes overexpressed in er- tumors  versus  <dig> of  <dig> in genes overexpressed in er+ tumors.
though the e2f site did not pass our multiple comparisons correction, published data support a role for these e2f sites in carcinogenesis.
prior efforts to identify a conditional regulatory program responsible for the coordinate regulation of sets of genes in multiple cancer types identified e2f as the lone tf universally overexpressed in multiple tumor types  <cit> .
the presence of e2f sites exclusively in genes overexpressed in er- brca tumors suggests that e2f plays a major role in this tumor type and may activate some target genes involved in cell cycle control  <cit> .
a caveat to our analyses is the realization that in some cases the motif count alone may not be considered to be a good predictor due to positional bias of a given motif relative to the transcriptional start site .
for some tfs, positional bias is likely to play a role in function.
for example, the motif tataaatw , well known for interactions with the basal transcription apparatus, shows a strong bias  <dig> bp upstream of the tss.
this spatial restriction is likely due to necessary interactions with the basal transcriptional apparatus   <cit> .
thus, motif copies present around - <dig> are likely to be functional while motifs distributed at other positions throughout the  <dig> kb upstream region would be predicted to be non-functional.
of our  <dig> phylogenetic motifs, only 32%  show positional bias, the majority of which are located within  <dig> bp of the tss.
the absence of any position bias for the vast majority of motifs in genes demonstrating disparate motif frequencies suggests a possible position-independent role in contributing to the observed expression patterns.
the lone phylogenetic motif showing significance, cagnygknaaa does not demonstrate positional bias.
a difficulty with any meta-analysis is that of study heterogeneity when one combines studies  <cit> .
meta-analyses on gene expression data are not immune from this criticism.
there are many factors that influence a designation of er+ and er- status in breast tumors, including assay sensitivity and the scoring system used.
the specific methods and assays for determining er+ and er- status are not available from oncomine and we were unable to account for this factor in our results.
many have proposed statistical methods for quantifying the heterogeneity in a meta-analysis data set  <cit> .
since heterogeneity manifests in an inflation of inter-study variance, a meta-analysis with any degree of heterogeneity tends to bias the effect size toward the null hypothesis  <cit>  and hence be conservative.
CONCLUSIONS
our meta-analysis was designed to identify genes showing consistent differences in gene expression patterns between er+ versus er- breast tumors.
the target genes identified provide a unified set of genes obtained across multiple analyses and their expression patterns may reflect the true biological complexity of breast tumors.
a small 10-gene meta-analysis signature to predict er status has recently been described  <cit> .
three genes identified in their study  overlap with our top 1% er+ upregulated genes.
from our results, a more highly refined set of gene targets can potentially be explored that would prove useful in the development of an improved biomarker assays for determining not only er status but also prognosis.
importantly, the overlap of  <dig> genes from our top 5% er- upregulated tumors with a set of  <dig> genes demonstrating overexpression in more than  <dig> types of undifferentiated cancers via meta-profiling identifies genes universally activated in cancer.
this list includes genes shown to be involved in the undifferentiated phenotype.
they include the melk kinase involved in mammalian embryogenesis, the apoptosis inhibitor birc <dig>  and multiple genes implicated in cell cycle control .
by screening the proximal promoter and 3'utr domains of our gene sets we wanted to identify both known tfbs, phylogenetically conserved motifs, and mirna seed sequences that differ in prevalence between er+ upregulated versus er- upregulated genes.
for any given site the disproportionate distribution between these gene sets may identify elements responsible for the co-regulation of groups of genes, and our analyses identified several significant elements in both the promoter and 3'utr regions.
moreover, er- genes had significantly shorter 3'utrs than er+ genes.
short 3'utrs are common for mirna anti-targets, which suggest that different mechanisms regulate groups of er+ and er- genes; that is, er+ genes may be mirna targets whereas er- genes may be anti-targets.
consistent with this hypothesis, er+ genes have significantly more putative mirna target sites in common than expected by 3'utr length alone, whereas er- genes have significantly less putative mirna target sites in common than expected by 3'utr length alone.
anti-target genes are commonly involved in basic cellular processes  <cit>  and in agreement with this, genes involved in the cell-cycle are significantly overrepresented in the er- genes .
clearly, our analysis is a starting point.
an examination of larger sequence domains upstream or these target genes may suggest additional elements showing differences in target abundance between these gene sets.
while our phylogenetic motifs were for the most part small , larger sequence elements such as enhancers that function at extended distances from these genes are likely to also play a role in the observed expression patterns.
the potential importance of promoter motifs in gene expression and of common polymorphisms that reside within these sites was highlighted by a recent survey of the promoter regions of nearly  <dig> genes in which 75% of the snps identified modify  putative tfbs  <cit> .
a survey of known polymorphisms  from existing databases  that reside within these motifs would also suggest the importance of these elements.
it would be of keen interest to explore if regulatory modules exist within these gene sets consisting of combinations of both known and phylogenetically conserved motifs.
approaches such as this have been described computationally for yeast, fly, mouse and humans  <cit> .
the recent use of comparative genomics tools from mammalian as well as evolutionarily distant species such as pufferfish  to identify phylogenetically conserved enhancers may also enable the identification of additional sequence elements responsible for the coordinate expression patterns seen for some of our genes  <cit> .
efforts such as these in conjunction with genomewide chromatin immunoprecipitation  studies of promoter regions will provide a more comprehensive view of the key elements modulating the observed gene expression patterns.
likewise, in the 3'utr, genomewide efforts to map snps to mirna target sites have revealed that many polymorphisms can either create new mirna target sites or can lead to their loss  <cit> .
genome-wide searches in humans have identified cis polymorphisms in putative mirna target sites that are likely contributors to phenotypic variation in humans and may to play a role in disease pathogenesis  <cit> .
future analyses will reveal whether snps in phylogenetically conserved promoter and 3'utr elements can influence breast cancer risk at the level of rna transcription or stability.

 variation in gene expression provides a quantifiable trait that has been employed to classify breast tumors. however it has long been known that the gene sets identified from independent laboratories fail to provide a unified set of genes thereby casting doubt on the biological implications of these profiles. despite these differences, two prognostic tests have recently been approved in the united states for clinical management of disease. from a diagnostic perspective, developing a unified gene profile that predicts both risk of recurrence and therapeutic response in diverse disease subtypes would be clinically useful. these gene sets could also provide an understanding of the mechanistic basis of malignancy. meta-analysis has been used as a formal summarization method in the clinical cancer literature for many years. recently, some groups have applied meta-analysis to gene expression microarrays. meta-analysis refers to a broad class of models used for summarizing and synthesizing studies to estimate their overall effect. rhoades, et al was among the first to demonstrate the usefulness of meta-analytic procedures on microarray data in prostate cancer. since then, there have been many contributions to the oncology literature by applying meta-analysis to microarrays, including breast cancer. one of the central goals in gene expression experiments is to identify the common regulatory themes and cis-elements responsible for the observed patterns of gene expression. this has been most successfully performed for the yeast saccharomyces cerevisiae where new regulatory genes have been suggested. however, metazoan expression patterns tend to be more complicated. one approach has been to combine expression data of orthologous genes from diverse organisms to build co-expression networks. in drosophilia gene networks have been proposed based upon the co-localization of tfbs with cis-regulatory modules. the availability of both mammalian and lower metazoan complete genome assemblies affords one the opportunity to identify phylogenetically conserved motifs in the array candidates. in addition to known tfbs, these phylogenetic motifs may identify important new cis-acting signals that modulate transcription or transcript stability and may be key elements in the observed expression patterns. a systematic comparison of both known and phylogenetic cis-elements between two sets of differentially expressed genes can serve to implicate these elements as common modulators in the observed gene expression patterns. our method incorporates a meta-analysis model to rank genes into groups of over- and under-expressed gene sets, based upon their relative importance between independent array studies. our analyses of gene expression patterns in er+ and er- breast tumors were performed across different array platforms on a diverse spectrum of patients. the two sets of genes showing the most disparate expression patterns between er+ and er- tumors provided an entry point with which to explore the possibility that specific sequence elements may be disproportionately represented in these two groups. we utilized known motifs in conjunction with comparative genomic resources to search for enriched dna elements in both the proximal promoter and 3'utr regions of these genes. our findings suggest that the differential gene expression in er+ vs. er- tumors may, in some cases, be mediated by specific sequence elements in either the promoter or 3'utr intervals. the motif distribution profiles between our gene sets identified both known and phylogenetically conserved elements that may play a role in these genes' co-expression.  forty-six percent of unique probes among the studies mapped many-to-one to unique unigene ids. the mean and median numbers of probes per unigene ids were and respectively. when we merged the studies in table for the meta-analysis data set, we retained the expression values for all probe combinations in all studies and this resulted in a multiplicative set of records in the database. approximately 12% of the unique ests in the oncomine database did not correspond to a unique unigene id. these were dropped from the analysis data sets. dfs, disease free survival; pr, progesterone status; ln, lymph node status; os, overall survival we focused our subsequent analyses on a select set of genes by taking medians across each unigene id's s/sd statistics. a scatter plot of the s/n versus abs/sd statistics appears in figure the distribution of the s/n values was bell-shaped with heavy tails. our criteria for selecting genes were to take the most extreme 1% and 5% values in both tails. we found it instructive to consider the ratio s/sd on the y-axis of figure where sd is the standard deviation of the addends of s. large values of this ratio indicate those genes with consistently significant p-values across all of the studies that we considered. the number of unigene ids with s/n scores in the top 1% and 5% were and respectively. the mean numbers of studies for genes present in the top 1% and top 5% classes were and respectively. our choice of reporting both top 1% and the top 5% for further screening was for crude management of false positives from bias correlated with each gene's relative ranking. many of the genes present in our top 1% upregulated list identified in our meta-analyses have previously been identified as overexpressed in er+ breast tumors, most notably the two transcription factors esr and gata our gene lists appear in additional file we next compared our top 1% and 5% upregulated gene lists in er+ and er- tumors to those prognostic genes utilized in the 70-gene signature associated with mammaprint® along with the 16-gene signature with the rt-pcr based oncotype dx® tests. although the array data defining the 70-gene profile was one of input datasets for our meta-analysis and the 16-gene signature datasets utilized two expression datasets were also employed for the analysis, we did not observe complete overlap in the genes identified. for the 70-gene signature our top 1% dataset identified an overlap of one and four genes respectively that were upregulated in er+ tumors versus er- tumors. only and genes overlapped in the top 5% dataset, respectively. alternatively, for the 16-gene signature, one and two genes, respectfully, from the top 1% gene sets were overexpressed in er+ versus er- tumors from our meta-analysis, while and genes, respectfully, overlapped in the top 5% list. differences in probes, arrays, and studies used in the meta-analyses may explain some of the differences between our gene lists and the gene lists from the two diagnostic tools. additionally, we compared our gene lists to a previously identified universal profile that uses genes overexpressed in a diverse spectrum undifferentiated cancers to predict neoplastic transformation. strikingly we observed only genes overexpressed from er- tumors to overlap with this gene signature. four genes overlapped in our top 1% list while genes overlapped from our top 5% list. these findings are in accord with the observation that er- tumors are more highly undifferentiated than er+ tumors and demonstrate more metastatic potential clinically. collectively these genes may identify a set of candidate genes predictive of metastatic potential in er- breast tumors. ingenuity pathway analyses we considered the relationship of our top 1% genes in the er+ and er- groups using ingenuity pathway analyses. our objective in using ingenuity was to characterize the functional role of our selected genes. ipa isolated genes for which it had documented associations, and created a series of networks based on the published literature. we were able to map of the genes comprising the sum of the 1% upregulated and 1% downregulated gene sets. from these networks, ipa queried its database of biological functions and scored each gene cluster with a p-value calculation. table shows the most common functions found among our most differentially-expressed genes. notably our top 1% genes upregulated in er- tumors contained genes showing association to cancer whereas only of the genes upregulated in er+ tumors were cancer-associated. promoter motif comparisons in dysregulated genes we tested the hypothesis that there was a significant difference in the occurrence of each motif between our two classes of genes using a fisher's exact test. we adjusted for multiple testing by applying the benjamini-hochberg p-value correction. we counted the number of genes in each class which were overexpressed in er+ tumors and contained a copy of each phylogenetic motif, and compared those to the number of genes overexpressed in er- tumors. for genes harboring multiple copies of a motif we counted these elements as a single motif event. we independently performed tests for both the top 1% and 5% of our genes. our initial query sets consisted of condensed transfac motifs and a second analysis comprised phylogenetically conserved mammalian promoter motifs as previously defined. sixty-nine of the phylogenetic motifs map to known tfbs defined in the transfac db v while represented novel phylogenetically conserved elements. we first examined whether any of known tfbs were disproportionately represented in our er+ and er- gene sets. abbreviated results appear in table while numerous motifs showed significance by fisher's exact testing only survived multiple testing correction. the first motif ktwgttt, a binding site for the sry transcription factor, was over-represented in the top 5% of er+ upregulated genes in the noncoding strand. for er+ overexpressed genes of genes contained the site while of er- overexpressed genes contained the site. the second site, abwcaggtrnr, a binding site for areb, was over-represented in the top 1% of er+ upregulated genes when both coding and noncoding strands were surveyed and contains an embedded e-box motif. twenty-five genes bore tcf sites in either strand amongst er+ upregulated genes while only genes contained the site amongst er- upregulated genes. the presence of tcf sites in nearly four times as many er+ upregulated genes versus er- upregulated genes may be an indirect mechanism for gene activation in er+ breast tumors. tcf has been shown to be induced by estrogen which in turn activates a cascade of downstream genes. additionally, the transcriptional repression of e-cadherin by tcf has been shown to lead to loss of the epithelial phenotype suggesting a role for this tf in late-stage carcinogenesis. we note that although e-cadherin was not identified in our meta-analysis, related genes, cdh and pcdh both of which lie in the top 5% of er- overexpressed genes, may be responsive to repression by tcf the over-representation of tcf binding sites in both strands of our top 1% genes er+ overexpressed tumors suggests that tcf may act as a transcriptional activator for these genes yet act as a transcriptional repressor in er- overexpressed genes. in addition to known sites, we sought to identify potential new regulatory motifs by examining the coding and noncoding strands with previously identified phylogenetic motifs in the top 1% and 5% of our s+ vs. s- genes. eleven of these motifs represented palindromic sequences and were scanned in only the coding strand when both strands were analyzed. again, while numerous motifs showed significance by fisher's exact testing only survived multiple testing correction. abbreviated results appear in table a single motif showed a significant difference between the er+ upregulated genes versus the er- upregulated genes when the non-coding strand was examined in our top 1% gene list. nineteen of er+ overexpressed genes contained at least copy of the motif while only of genes contained the motif in the er- overexpressed genes. this phylogenetic motif does not map to any known tfbs and represents a new target for exploration. *center of motif relative to transcriptional start site analysis of the 3'utr we next screened for regulatory elements in the 3'utr of our genes sets. less is known about functional motifs in 3'utrs than about functional motifs in promoter regions, but evolutionary conserved motifs in 3'utrs may, as in promoter regions, indicate regulatory sites. we therefore used a previously identified set of evolutionary conserved 3'utr motifs. although the function of half of these motifs is unknown, the remaining half has a/t rich elements believed to be involved in mrna stability or represent likely micrornas binding sites. we used the same refseq id's to harvest the annotated 3'utrs of our gene sets as described in the methods. surprisingly, we observed a significant difference in the median 3'utr lengths between our gene sets. the top 1% genes overexpressed in er+ tumors contained a median 3'utr length of kb, while genes overexpressed in er- tumors contained a median 3'utr length of kb. a similar trend was observed when we examined the top 5% of genes sets. in this set, er+ upregulated genes had a median utr length of kb while the er- genes had a length of kb. microrna target genes have longer 3'utrs, whereas anti-targets have shorter 3'utrs. thus, the difference in 3'utr length suggests a difference in mirna targeting prevalence between the er+ and er- genes. the most significant evolutionary conserved motif in the top 1% and top 5% genes correspond to a potential mirna target site; yactgccr and wgcctta have seed complementarity to mir-34/mir- and mir- the mirna seed region – nucleotides 2– from the 5' end – is the most important factor for mirna target site recognition. fisher's exact tests on the mirna seed site occurrence counts, corrected for multiple testing, seemingly confirm that the er+ genes are preferentially regulated by mirnas, as all the significant seeds are overrepresented in the er+ upregulated genes. there is, however, a potential problem with using the fisher's exact test for the 3'utr sets. if motif occurrences were random, we would expect the er+ genes to have more motif occurrences than the er- genes have, as the er+ genes have longer 3'utrs. thus, to determine whether there is a significant difference in mirna regulation between the er+ and er- genes, we had to address whether the occurrences of mirna seed sites in the two sets were significantly different from what we would expect by chance. we therefore ran a set of randomization experiments where we compared the observed number of seed site occurrences in the er+ and er- genes' 3'utrs with those in random gene sets that had 3'utr lengths similar to the er+ and er- 3'utrs. we found that all of the seeds identified by significant fisher's exact tests do occur significantly more frequently in the er+ 3'utrs than in 3'utrs from random gene sets. moreover, these seeds also occur significantly less frequently in the er- 3'utrs than in random gene sets. thus, it seems that whereas several mirnas may coordinately regulate some of the er+ genes, some of the er- genes may collectively avoid being regulated by the same mirnas. previous studies have identified several mirnas that are aberrantly expressed in breast cancers. together the aberrantly expressed mirnas in these studies mapped to unique mer seed sequences of which three were among our ten most significant mer motifs. the three corresponding mirnas are all overexpressed in breast cancers. none of the ten most significant mer motifs are from mirnas reported to be differentially expressed in er+ and er- tumors; the most significant mer is ranked 25th, is significantly more abundant in er+ genes than expected by random, and is from mir- which is downregulated in er+ tumors. discussion micrornas are small noncoding rnas that recognize complementary target sequences in mrnas and prompt either translational repression or rna degradation. micrornas play important roles in cancer. iorio et al., for example, recently revealed that deregulation of multiple mirnas can be correlated to pathogenic features such as estrogen or progesterone receptor status and tumor stage for breast cancers. in addition, shorter postoperative survival times for patients with lung tumors can be predicted by measuring mirna let-. thus mirnas can be used both as classifiers of breast tumor type and as predictors of survival of lung cancer patients. micrornas preferentially target 3'utrs that have short sequences with perfect complementarity to nucleotides 2– or 2– in the mirna's 5' region – the seed region. as mirna regulation may explain gene co-expression, we therefore included the mer and mer seed sequences for all human mirna sequences known at the time of the study. we note that not all known human mirnas are highly evolutionary conserved and these seed sequences therefore supplement the mirna-related evolutionary conserved motifs. since we identified sets of genes that demonstrated differential expression between er+ and er- tumors, we reasoned that some of these genes may contain common cis-regulatory motifs contributing to their co-regulation. we would predict that these sites may, in some cases, be disproportionately represented between genes upregulated in er+ tumors versus genes upregulated in er- tumors perhaps allowing one to identify genes sharing common regulatory pathway. computational tools exist to identify tfbs based upon over-representation of conserved motifs in datasets. other approaches aim to identify transcription factors which bind to tfbs based on the relatedness of expression profiles between the tf and the target genes they are postulated to regulate. a combined approach utilizing expression measurements of tissue-specific gene sets in conjunction with orthologous tfs from humans and mouse provides for enhanced accuracy in predicting bone fide cis-regulatory elements. for the most part these searches are guided by biologically confirmed tfbs interactions identified in the transfac database; however, this approach may fail to identify motifs that may be evolutionarily conserved amongst mammals. in addition to known sites that remained significant after multiple testing correction, many additional sites, and their associated transcription factors, warrant comment. a second important tfbs, ctttga, the binding site for lymphoid enhancer-binding factor, in the top 1% coding strand er+ overexpressed genes, failed rigorous multiple testing where of genes contained ≥ site versus of genes in er- gene set in table nonetheless there is strong biological evidence supporting the role of lef in tumorogenesis. the lef binding site ctttga is one of the primary binding sites in the wnt signaling pathway which regulates cell-cell adhesion and many morphogenetic events during mammary development and possibly cancer binding of wnt proteins with frizzled protein prevents degradation of β-catenin, which subsequently translocates to the nucleus and binds transcription factors of the tcf/lef family. several tumors are known to have an altered β-catenin signaling pathway including colorectal and lymphoblastic tumors. mutations in the wnt pathway genes can result in β-catenin stabilization and activation of lef/tcf-induced transcription. recent studies have demonstrated sebaceous tumors harboring lef mutations interfere with β-catenin-binding domain of lef and transcriptional activation. common human carcinomas also carry mutations in the β-catenin-binding domain of lef. our data suggest that mutations in lef or tcf binding sites in genes that inactivate wnt signaling could contribute to breast tumorogenesis. we did not find the estrogen receptor binding site over-enriched in any our analyses. this is not surprising as our survey was confined to the immediate kb promoter region. we point out that estrogen may be playing an indirect role on genes in er+ overexpressing tumors via the activation of tf such as tcf which in turn activate downstream targets. additionally, it is possible that differences in er binding sites do exist between our gene sets but these sites may reside at distances much further upstream. recent reports indicate that only two-thirds of er tfbs can be localized to the proximal promoter region of rna polymerase ii genes. we also note that the e2f binding site consistently ranked amongst the top motifs identified when screening the non-coding strand. in the non-coding strand of the top 1% gene sets, more e2f sites were observed in genes overexpressed in er- tumors versus of in genes overexpressed in er+ tumors. though the e2f site did not pass our multiple comparisons correction, published data support a role for these e2f sites in carcinogenesis. prior efforts to identify a conditional regulatory program responsible for the coordinate regulation of sets of genes in multiple cancer types identified e2f as the lone tf universally overexpressed in multiple tumor types. the presence of e2f sites exclusively in genes overexpressed in er- brca tumors suggests that e2f plays a major role in this tumor type and may activate some target genes involved in cell cycle control. a caveat to our analyses is the realization that in some cases the motif count alone may not be considered to be a good predictor due to positional bias of a given motif relative to the transcriptional start site. for some tfs, positional bias is likely to play a role in function. for example, the motif tataaatw, well known for interactions with the basal transcription apparatus, shows a strong bias bp upstream of the tss. this spatial restriction is likely due to necessary interactions with the basal transcriptional apparatus. thus, motif copies present around - are likely to be functional while motifs distributed at other positions throughout the kb upstream region would be predicted to be non-functional. of our phylogenetic motifs, only 32% show positional bias, the majority of which are located within bp of the tss. the absence of any position bias for the vast majority of motifs in genes demonstrating disparate motif frequencies suggests a possible position-independent role in contributing to the observed expression patterns. the lone phylogenetic motif showing significance, cagnygknaaa does not demonstrate positional bias. a difficulty with any meta-analysis is that of study heterogeneity when one combines studies. meta-analyses on gene expression data are not immune from this criticism. there are many factors that influence a designation of er+ and er- status in breast tumors, including assay sensitivity and the scoring system used. the specific methods and assays for determining er+ and er- status are not available from oncomine and we were unable to account for this factor in our results. many have proposed statistical methods for quantifying the heterogeneity in a meta-analysis data set. since heterogeneity manifests in an inflation of inter-study variance, a meta-analysis with any degree of heterogeneity tends to bias the effect size toward the null hypothesis and hence be conservative.  our meta-analysis was designed to identify genes showing consistent differences in gene expression patterns between er+ versus er- breast tumors. the target genes identified provide a unified set of genes obtained across multiple analyses and their expression patterns may reflect the true biological complexity of breast tumors. a small 10-gene meta-analysis signature to predict er status has recently been described. three genes identified in their study overlap with our top 1% er+ upregulated genes. from our results, a more highly refined set of gene targets can potentially be explored that would prove useful in the development of an improved biomarker assays for determining not only er status but also prognosis. importantly, the overlap of genes from our top 5% er- upregulated tumors with a set of genes demonstrating overexpression in more than types of undifferentiated cancers via meta-profiling identifies genes universally activated in cancer. this list includes genes shown to be involved in the undifferentiated phenotype. they include the melk kinase involved in mammalian embryogenesis, the apoptosis inhibitor birc and multiple genes implicated in cell cycle control. by screening the proximal promoter and 3'utr domains of our gene sets we wanted to identify both known tfbs, phylogenetically conserved motifs, and mirna seed sequences that differ in prevalence between er+ upregulated versus er- upregulated genes. for any given site the disproportionate distribution between these gene sets may identify elements responsible for the co-regulation of groups of genes, and our analyses identified several significant elements in both the promoter and 3'utr regions. moreover, er- genes had significantly shorter 3'utrs than er+ genes. short 3'utrs are common for mirna anti-targets, which suggest that different mechanisms regulate groups of er+ and er- genes; that is, er+ genes may be mirna targets whereas er- genes may be anti-targets. consistent with this hypothesis, er+ genes have significantly more putative mirna target sites in common than expected by 3'utr length alone, whereas er- genes have significantly less putative mirna target sites in common than expected by 3'utr length alone. anti-target genes are commonly involved in basic cellular processes and in agreement with this, genes involved in the cell-cycle are significantly overrepresented in the er- genes. clearly, our analysis is a starting point. an examination of larger sequence domains upstream or these target genes may suggest additional elements showing differences in target abundance between these gene sets. while our phylogenetic motifs were for the most part small, larger sequence elements such as enhancers that function at extended distances from these genes are likely to also play a role in the observed expression patterns. the potential importance of promoter motifs in gene expression and of common polymorphisms that reside within these sites was highlighted by a recent survey of the promoter regions of nearly genes in which 75% of the snps identified modify putative tfbs. a survey of known polymorphisms from existing databases that reside within these motifs would also suggest the importance of these elements. it would be of keen interest to explore if regulatory modules exist within these gene sets consisting of combinations of both known and phylogenetically conserved motifs. approaches such as this have been described computationally for yeast, fly, mouse and humans. the recent use of comparative genomics tools from mammalian as well as evolutionarily distant species such as pufferfish to identify phylogenetically conserved enhancers may also enable the identification of additional sequence elements responsible for the coordinate expression patterns seen for some of our genes. efforts such as these in conjunction with genomewide chromatin immunoprecipitation studies of promoter regions will provide a more comprehensive view of the key elements modulating the observed gene expression patterns. likewise, in the 3'utr, genomewide efforts to map snps to mirna target sites have revealed that many polymorphisms can either create new mirna target sites or can lead to their loss. genome-wide searches in humans have identified cis polymorphisms in putative mirna target sites that are likely contributors to phenotypic variation in humans and may to play a role in disease pathogenesis. future analyses will reveal whether snps in phylogenetically conserved promoter and 3'utr elements can influence breast cancer risk at the level of rna transcription or stability. 
we surveyed both the proximal promoter and 3' untranslated regions of our top-ranking genes in each expression group to test whether common sequence elements may contribute to the observed expression patterns. many of the genes identified as overexpressed in er- tumors were previously identified as expression markers for neoplastic transformation in multiple human cancers. this allows one to compare expression data from a diverse collection of breast tumor array datasets generated on either cdna or oligonucleotide arrays. we gathered expression data from published microarray studies examining estrogen receptor positive and estrogen receptor negative brca tumor cases from the oncomine database. utilizing a combination of known transcription factor binding sites, evolutionarily conserved mammalian promoter and 3'utr motifs, and microrna seed sequences, we identified numerous motifs that were disproportionately represented between the two gene classes suggesting a common regulatory network for the observed gene expression patterns. moreover, our motif analysis identified a collection of specific cis-acting target sites which may collectively play a role in the differential gene expression patterns observed in er+ versus er- breast cancer tumors. importantly, the gene sets and associated dna motifs provide a starting point with which to explore the mechanistic basis for the observed expression patterns in breast tumors.
 gene expression measurements from breast cancer tumors are established clinical predictive tools to identify tumor subtypes, identify patients showing poor/good prognosis, and identify patients likely to have disease recurrence. however, diverse breast cancer datasets in conjunction with diagnostic clinical arrays show little overlap in the sets of genes identified. one approach to identify a set of consistently dysregulated candidate genes in these tumors is to employ meta-analysis of multiple independent microarray datasets. this allows one to compare expression data from a diverse collection of breast tumor array datasets generated on either cdna or oligonucleotide arrays.  we gathered expression data from published microarray studies examining estrogen receptor positive and estrogen receptor negative brca tumor cases from the oncomine database. we performed a meta-analysis and identified genes that were universally up or down regulated with respect to er+ versus er- tumor status. we surveyed both the proximal promoter and 3' untranslated regions of our top-ranking genes in each expression group to test whether common sequence elements may contribute to the observed expression patterns. utilizing a combination of known transcription factor binding sites, evolutionarily conserved mammalian promoter and 3'utr motifs, and microrna seed sequences, we identified numerous motifs that were disproportionately represented between the two gene classes suggesting a common regulatory network for the observed gene expression patterns.  some of the genes we identified distinguish key transcripts previously seen in array studies, while others are newly defined. many of the genes identified as overexpressed in er- tumors were previously identified as expression markers for neoplastic transformation in multiple human cancers. moreover, our motif analysis identified a collection of specific cis-acting target sites which may collectively play a role in the differential gene expression patterns observed in er+ versus er- breast cancer tumors. importantly, the gene sets and associated dna motifs provide a starting point with which to explore the mechanistic basis for the observed expression patterns in breast tumors. 