A comparative genomics approach to identifying the plasticity transcriptome

Developing binding site consensus sequences

Previous analyses of transcription factor binding sites have suffered from low degeneracy based upon comparison of a target sequence to a single overrepresented binding site or high degeneracy and a high false positive rate based upon inclusion of related IEG subfamily members with similar but distinct sequence specificities. When related transcription factors (such as zif268, egr-2, and egr-3 or c-fos/c-jun and FosB/JunD heterodimers) with overlapping binding site consensus sequences are pooled, binding site consensus is significantly relaxed and false positive rates are increased. Our method sought to reduce the rate of false positives by using a smaller number of experimentally verified, high-quality transcription factor binding sites using frequency matrices that have been experimentally developed for zif268 [45, 48] and CREB [44, 48]. A schematic of the consensus sequences used is shown in Fig. 1.

An analogous search was attempted with AP-1, whose components, fos and jun, are also IEGs. The resulting predicted targets were not enriched for promoter regions and showed little conservation across species, facts we attribute to excessive degeneracy of the known AP-1 binding sequences used in training the computational methods. While we believe this negative result helps to reinforce the significance of the positive CREB and zif268 results, we have omitted further discussion of AP-1 in this report.

Computational genomic analysis: mouse, rat, and human genomes

In order to identify conserved gene candidates subject to activity-dependent regulation, the entire database of mouse, rat, and human promoter sequences was searched for CREB and zif268 binding sites. In total, 18,071 mouse genes, 5,943 rat genes, and 19,794 human genes were examined (summarized in Table 1; full gene lists can be found in Additional Files 1, 2, 3, 4). This represents all the annotated promoters from these species (approximately 50% of the estimated distinct mouse coding regions and two-thirds of the total estimated distinct human genes; [49–51]). CREB binding sites were predicted in 6% of mouse promoters, 7% of human promoters, and 11% of rat promoters (Fig. 2a). Zif268 binding sites were predicted in 8% of mouse promoters, 6% of human promoters and 4.6% of rat promoters (Fig. 2b). The precise sequence and location of binding site for each species is annotated in Additional Files 5, 6, 7. Because the mouse and human gene sets were more complete and representative of the total number of coding regions within those species than found in the rat dataset, we will primarily refer to the human and mouse searches from this point onward.

As we were interested in using comparative genomics to improve hit quality, we also examined those genes with homologues in the mouse and human datasets, a total of 13,365 genes (homologene dataset; [52]). In general, the frequencies of binding sites predicted for genes in the homologene dataset for either species in isolation were similar to those found for non-homologous genes in the same genomes (Table 1), suggesting that the homologous genes are a representative set of genes with respect to binding site predictions. We then examined the frequencies of binding site predictions conserved across both members of a gene pair. The number of hits from this search was lower than observed for individual genomes: in the set of mouse-human homologues, 2.66% have a conserved CREB site and 1.24% have a conserved zif268 site (Fig. 2c).

Estimating rate of false positives

To develop a rigorous estimate of the quality of identified targets, we applied two methods for calculating their positive predictive value, defined as the probability a predicted hit is correct (see Methods for a complete description). This measure provides a conservative estimate of the percentage of identified transcription factor binding sites that we expect to be functional binding sites. For this purpose, we examined binding site location and cross-species conservation.

Fig. 3 shows an analysis of position specificity of binding site predictions within the promoter regions. A pronounced peak in the relative frequency of binding site locations was observed for both CREB (Fig. 3a; see also [46]) and zif268 (Fig. 3b). Both binding sites were more than five times more likely to occur in the 50 bp closest to the annotated start site compared to more distant or intergenic sequence. This frequency increase was observed in analysis of the mouse and human gene datasets, as well as in the conserved mouse and conserved human gene datasets (Fig. 3c and 3d). Furthermore, the low frequency of binding sites in regions most distal to the annotated start site (i.e. 800 bp upstream) was comparable to the frequency in intergenic regions (see Methods), suggesting that our selection of a 1,000 bp region surrounding transcription start was sufficient to capture most of the meaningful binding sites that might regulate gene expression. Using intergenic (i.e. 50,000 nt upstream of the annotated start site) hit frequencies to estimate false positive rate, we derived a positive predictive value of 0.64 for CREB hits in human and mouse promoters. Zif268 had a slightly higher positive predictive value of 0.81 in mouse and 0.74 in human (Table 1). The full set of predicted binding site locations is provided in Additional Files 1, 2, 3, which present the gene, location relative to transcription start, and sequence of each putative binding site for the mouse, human, and rat genomes, respectively.

Comparative genomics was used to provide a second method of estimating the positive predictive value. In the set of mouse-human homologues, 2.7% have a conserved CREB site and 1.2% have a conserved zif268 site. We proposed that the set of homologue pairs with conserved binding sites would provide a conservative set of true positive hits. The positive predictive value of the hits could then be estimated from the excess of conserved hits beyond what would be predicted from the hit rates in the individual genomes on the assumption of independence between genomes (see Methods). In contrast to the results using location specificity, the conservation approach yielded higher estimates of positive predictive value for CREB (0.83) than for zif268 (0.56). This difference is most likely to due changes in zif268 targets or promoter GC content between species. We can further validate these results by measuring conservation in binding site positions for the conserved hits. Comparing the positions of the conserved predicted sites between mouse and human yields correlation coefficients of 0.18 for CREB (p-value < 0.001) and 0.21 for zif268 (p-value < 0.01). There is considerable noise in location specificity, which might be explained by insertions or deletions in the promoters since the separation of human and mouse lineages or by errors in annotated transcription start sites in either species. The positions nonetheless show significantly greater correspondence than can be explained by chance.

Overlap of CREB and zif268 target genes

Prior analyses of activity-dependent gene expression have not addressed whether the transcription factors selected for analysis regulate separate or partially overlapping target genes. We reasoned that genes containing binding sites for both transcription factors might be more likely to have an essential role in mediating activity-dependent changes in gene expression. To examine this, we looked at the relative frequencies of genes with both CREB and zif268 binding sites in the mouse, human, and rat gene datasets compared to the frequencies of genes with either site in isolation. Target genes were generally non-overlapping, suggesting that each transcription factor may regulate a distinct subset of genes, a possibility that might allow a cell greater combinatorial control over stimulus-specific transcriptional programs.

Intriguingly, however, the amount of overlap between zif268 and CREB targets was greater than would be expected by chance, suggesting that at least a subset of target genes may be coregulated by the two transcription factors analyzed. Statistical analysis showed that the predicted binding sites shared more targets in common than would be expected by chance for both mouse and human genomes, with the results weakly significant for mouse (p < 0.05), but not significant for human (p = 0.18).

Only six genes with conserved binding sites for both transcription factors were found. A significant fraction of these common targets were transcription factors. FosB (Fig. 4a), Jund1 (Fig. 4b), and Maff (Fig. 4c) are all members of the AP-1 family of transcription factors. The Skil (Fig. 4d) transcription factor is a member of the SKI/SNO/DAC family which are known to associate with AP-1 under some conditions [53]. The observation that a specific group of transcription factors can be regulated by both CREB and zif268 implicate these genes in transcriptional networks of activity-regulated gene expression.

Computational analysis of target gene sets

We sought to characterize the functional properties of the derived gene set through an unbiased computational search for functional gene classes significantly over- or under-represented in our hit set. We conducted this analysis by applying the GOstat web resource [54] to the set of CREB and zif268 conserved and species-specific targets. Because there were relatively few conserved targets in the homologene datset, GOstat analysis found few hits, such as RNA processing and localization, only for CREB target genes.

Because of the increased statistical power of using more data (versus higher quality data), we chose to further this analysis using the species-specific gene target lists. We chose to present the mouse target list because the comparative data available has primarily been carried out in rodents. As with the conserved gene list, genes with CREB consensus sites showed significant overrepresentation for targets involved in RNA processing, but were underrepresented in electrophysiologically important transmembrane receptor targets (Additional File 8). However, the few receptor and channel targets identified in this analysis may have critical functional importance. For example, seizure-dependent changes in expression of three targets, Kcnk1, Kcnmb4 and Hcn2, have been found [55], and these two genes represent targets for neuroprotective or anticonvulsant agents. Targets of zif268 are also underrepresented for transmembrane receptors, but are overrepresented for transcription factors as well as genes with neural-specific functions such as neuron development and axonogenesis. A full list of these targets can be found in Additional File 9.

Comparison to experimental analyses

Although there are few well-documented differences between the transcriptional regulation of homologous genes in mouse and humans, our analysis was successful at identifying one of the few well-characterized mutations in a CREB binding site between a mouse and human gene that influence transcription, the glycoprotein hormone alpha subunit [62]. This gene carries an identified CREB site in human, and a single gene mutation in the CREB binding site in mouse abolishes placental expression of this gene. Accordingly, this gene was identified as a human but not a mouse CREB target gene.

Several experimental studies have sought to characterize CREB or zif268 targets by searching for genes upregulated after transcription factor activation. We compare the genes isolated in a subset of these studies to the targets identified here. Microarray results of regulated genes after 1–5 weeks of hippocampal overexpression of a constitutively active form of CREB [37] were compared to the results from our search, using a program that takes into account relative changes in transcript levels [63]. We found that predicted CREB targets from our search were overrepresented in that microarray dataset, but by a statistically insignificant amount. An analysis of CREB targets from a similar computational search [64]showed slightly less enrichment in this microarray dataset, but this enrichment was also not significant. The lack of significance in overlap is likely due to the extended duration of VP-16 CREB overexpression, which complicates an interpretation of this result, as well as the fact that we searched specifically for direct targets of CREB while a microarray study would be expected to capture both direct and downstream targets.

We also compared overlap between zif268 targets identified in our computational analysis and a recent zif268 overexpression study [49]. We observed an extremely low rate of overlap between these two studies, where only 6.7% (9/135) regulated genes were represented in either the mouse or human zif268 target lists generated in our study (Table 3). Only one gene, the small glutamine-rich tetratricopeptide, showed a conserved zif268 site in both mouse and human as well as regulated mRNA levels after zif268 overexpression. The relatively low level of target overlap in this and other studies may be due to the duration of transcription factor overexpression (days) and ensuing cascades of altered gene expression (i.e., not all targets are directly regulated by the overexpressed transcription factor). Our decision to favor a high-specificity gene list at the cost of lower sensitivity was likely also a factor. While the exact criterion used by James et al. study [57] to define target sites was unclear, likely differences in the target specificity/sensitivity trade-off would be expected to yield poor correspondence between the data sets.

Comparison to other computational analysis

A similar comparative genomics analysis was carried out by Conkright et al., using a hidden Markov search algorithm trained on ten well-characterized CREB binding sites [64]. Of 78 conserved CREB targets identified by Conkright et al. and 356 identified by our study, 25 were common to both studies (Table 2). Although this is a strong overlap relative to what one would expect by chance, it is nonetheless curious that it was not higher. We attribute this fact to likely assignment errors in both sets as well as likely differences in annotation. We believe our study is likely to have yielded a higher quality set of predicted binding sites based on the fact that we have access to more recent genome annotations, search in a more tightly focused region (1.2 kb versus 10 kb), search relative to transcription rather than translation start, and use a prediction algorithm that would screen out some possible spurious predictions likely with a hidden Markov model approach. This methodological argument that our gene set is closer to the ground truth is supported by the fact that while the two studies predict comparable numbers of CREB sites in mouse and human individually (1050 mouse and 1389 human for the present study and 1349 mouse and 1663 human for Conkright et al.) our predicted sites were validated by cross-species conservation at a rate several-fold higher (356 validated versus 78). The similarities of our gene set to that of Conkright et al. thus provide good validation that both approaches find meaningful gene sets, but the deviations do not challenge the accuracy of our set.

Promoter database

We compiled a database of gene promoter regions using sequences from mouse build mm6 [71], rat build rn3 [72], and human build hg17 [49, 51] of the UCSC Genome Bioinformatics Resource [73]. Transcription start sites for these promoters organized by mRNA accession number were found in the table "knownGene.txt" for each build. Where promoter regions were reported within 50 bp of each other, only the one earlier on the chromosome was used, as the copies were presumed to be duplicates of the same promoter region (derived from otherwise identical mRNAs of different lengths). Incomplete promoters with missing sequence data were also removed from the analysis. Annotated promoters included both TATA-box containing and TATA-less genes. The full promoter list was annotated with gene name, symbol, and accession number using the NCBI gene resources [52, 74]. In total, 18,071 mouse promoters, 19,794 human promoters, and 5,943 rat promoters were analyzed (Tables 1, Additional files 1, 2, 3).

When searching for candidate genes, we defined a putative promoter to be the genetic sequence from -1,000 bp to +200 bp of each transcription start for human, mouse, and rat genes. A set of intergenic sequences was also compiled for human and mouse to construct a "random" control dataset of 1,200 bp sequences using the regions from -51,200 bp to -50,000 bp relative to each transcription start site, where the transcription factors CREB and zif268 are not likely to have regulatory function. Due to a decrease in sequence quality further away from transcription start, distal sequence regions were available for only 77% of total genes, leaving 13,475 mouse intergenic regions and 15,178 human intergenic regions far analysis (Table 1). In order to confirm location-specificity trends inferred from the 1,200 bp regions, an additional search was run for each gene on an extended promoter region (-6,000 bp to +200 bp). We saw no significant difference in the region between -1,000 bp and -6,000 bp compared to the -51,200 bp to -50,000 bp region, suggesting that the initial search had identified the majority of sites with likely function.

The Homologene database provided us with human and mouse homologous pairs based on gene accession number [52], yielding 13,365 homologous gene pairs (Table 1, Additional File 4). A binding site prediction was defined as conserved if the same binding site type was predicted in the promoters of both homologous genes, without regard for position in the promoter.

Transcription factor binding site inference

The goal of this search was to identify a transcription factor binding site compared to its background. This is the case when the probability that it is a binding site is greater than the probability that the sequence would be observed by chance.

P(Model) > P(Background)

log(P(Model)) > log(P(Background))

log(P(Model)) - log(P(Background)) > 0

Position specific scoring matrices (PSSMs) or position weight matrices (PWMs), are a well-established method of motif finding [75, 76]. We used a variant of them to find the log probability of a sequence being a part of the model. These methods are similar to the transcription factor binding site search available through the database of transcription start sites [77]. Binding site frequency matrices for CREB and zif268 were obtained from the Transfac [48, 78] public database (see Fig. 1). These matrices give the frequency of each nucleotide in each position of the binding site. Scoring matrices for the present study were created from the Transfac frequency matrices with the following equation:

$S_{n,p}=\log (\frac{A_{n,p}+b}{{\displaystyle {\sum }_i{A}_{i,p}+4b}})$

where S is the scoring matrix, A is the frequency matrix, n is the nucleotide, p is the position within each binding site. The pseudocount, b, is set at the relatively small value of 0.25 to allow limited tolerance of base-pairs which have never been observed in a given position for a binding site. When comparing this scoring matrix to a sequence of the same size, adding the scores for the nucleotide that is at that same position in the sequence gives you the log of probability that the sequence matches the model.

The probability that a sequence is not a binding site is based on background dinucleotide frequencies. For each individual species, we went through all promoters and calculated the probability of each dinucleotide transition. For instance, $P(AT)=\frac{Observed(AT)}{Observed(Nucleotides)}$. We also calculated the probability of observing each nucleotide individually. The probability of observing any sequence can be calculated from those probabilities by multiplying the probability of the first nucleotide by the probability of each nucleotide transition. The log probability can be found by adding the log of each probability.

$\begin{array}{l}P(ACT)=P(A)*\frac{P(AC)}{P(A)}*\frac{P(CT)}{P(C)}\\ \log (P(ACT))=\log (P(A))+\log (\frac{P(AC)}{P(A)})+\log (\frac{P(CT)}{P(C)})\end{array}$

The goal of this study was to create a comprehensive list of possible transcription factor targets. The log of the sequence length is often subtracted to correct for the number of possible sites being searched. While subtracting by the full log of the sequence length would provide a more rigorous control, we deliberately chose to increase the sensitivity of the method at the expense of specificity. This increases the number of targets found while decreasing the average quality of the target genes, allowing us to better take advantage of homology to eliminate false positives. Instead of subtracting the log of the sequence length, we subtracted a smaller "correction value". The larger the correction value, the more stringent the search is. To analyze the effects of decreasing specificity on the quality of the target gene set, we use a measure called the positive predicted value. Intuitively, this measure is the probability that a predicted site is a true positive. The positive predictive value is defined as $\frac{TruePositives}{FalsePositives+TruePositives}$.

In terms of our data, the positive predictive value is

The expected number of sites is the number of sites expected to be conserved if there was no association between a binding site existing in mouse and its human homologue. A series of possible correction values are plotted against the positive predicted value (Additional File 10). Because it is the point at which CREB and zif268 positive predictive values plateau, we chose to use a correction value of 300. The final positive predictive value based on comparative genomics is found in Table 1 under "Homologues." A binding site is considered a hit if the final calculated score is above zero. The equation used to determine the final score is given below:

Score = log(P(Model)) - log(P(Background)) - log(correction_value)

Global data analysis

The positive predictive values for the individual species is calculated by comparing the promoter region to the intergenic region (see Table 1). It is still calculated as $\frac{ObservedSites-ExpectedSites}{ObservedSites}$, but the observed sites is now the percentage of promoter targets with a binding site and the expected sites is the percentage of intergenic regions with that binding site.

Associations between co-occurring binding sites were analyzed by applying a 2-tailed Fisher's exact test [79], using a web-based calculator [80], to the 2 by 2 contingency table of counts of occurrence of either, both, or neither site.

Analysis of the function of the activity-dependent transcription factor targets was done using GOstat [81], an online tool for finding overrepresented ontologies in a set of genes [54]. The list of targets for each transcription factor binding site and species was searched against the entire list of promoters for overrepresentation of different gene ontology classes.

The VP16-CREB expression data [37] was obtained from the Gene Expression Omnibus (Barrett), record GSE3965, on the NCBI website. The SOFT data files were converted to a single matrix using the GEOquery as part of the Bioconductor [82] package for R. We separated out the expression data into two groups: a control group where CREB is expressed at normal levels (On dox, rev, wt) and an experimental group where VP16-CREB has been active from one to five weeks (1w, 2w, 5w) [37]. Only experiments where the entire hippocampus is dissected are used, not micro-dissected CA1 regions of the hippocampus, where region-specific gene expression could bias the results. The data was analyzed using the GSEA software package for Windows [64]. The change in gene expression levels between the control and the super active CREB were determined by a signal-to-noise metric. We created genesets out of our conserved CREB targets and the conserved CREB targets identified by Conkright et al. [65]. Symbols present on these lists but with no corresponding microarray probe are ignored. The expression levels corresponding to the symbols in the genesets are queried for enrichment in the control versus the experimental microarray datasets using the methods described in Subramanian et al. A list of genes identified as the leading edge subset [64], which are genes that contribute to the enrichment of the CREB targets in the microarray data, are listed as the overlap between our dataset and VP16-CREB dataset [37].